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ABSTRACT 


The subject study addressed computational aspects of (l) flutter opti- 
mization (minimization of structural mass subject to specified flutter 
requirements), (2) methods for solving the flutter equation, and (3) efficient 
methods for computing generalized aerodynamic force coefficients in the 
repetitive analysis environment of conqDuter-aided structural design. The 
principal results of the study are summarized in a companion report. The pre- 
sent report contains supplemental data and supporting information on various 
aspects of the study including the following: 

Details of a two-dimensional Regula Falsi approach to solving the 
generalized flutter equation are presented. 

The method of Incremented Flutter Analysis and some of its applications 
are discussed. 

The use of Velocity Potential Influence Coefficients in a five-matrix 
product formulation of the generalized aerodynamic force coefficients is 
delineated. Options for computational operations required to generate gen- 
eralized aerodynamic force coefficients are compared in detail. 

Theoretical considerations related to optimization with one or more 
flutter constraints are presented as well as practical e::^erience with an 
actual structural design problem. 

Expressions for derivatives of flutter-related quantities with respect 
to design variables are presented. Included are flutter-speed derivatives 
with variable modes . 

A bibliography is included. 
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SYMBOLS AND DEFINITIONS 


[ ] 


{ } 

L J ^ 
C J 


square, rectangular matrix 

transpose of a matrix 
column matrix 
ro:.: /matrix 
diagonal matrix 
determinant 


[XI] 


matrix [x] with last column omitted 


[X] 


matrix [X] with last row omitted 


A 


constant multiplying the additive column of design 
variable increments (equation 5 * 17 ) 


a ,a ' amplitudes of successive cycles 

n n+1 


[A(ik)] 

[A(p)j 


aerodynamics matrix (function of k and Mach number), 
modalized aerodynamics matrix 

aerodynamics matrix (function of p and Mach number) 


[A(p)] 

[Aj] 


|pV^[A(p)] " 

j = 0,1, 2, 3; defines [A(k)] in cubic spline interpolation 
(equation ^.37) 



A.j0c) 


generalized aerodynamic force coefficients 


[AIC] ,[AIC(k)] basic aerodynamics influence coefficients (function of k 

and Mach number) defined by equation 4^1 | 

[AIC.] j = 0,1. 2.3; defines [AlC(k)] in cubic spline interpolation 

■ (similar to equation 4.37) 

[AW ] j = 0,1, 2, 3; defines [AW(k)] = [AlC('k) ][W(k) ] in cubic 

spline interpolation (similar to equation 4.37) 

B constant multiplying the subtractive column of design 

variable increments (equation 5«I7) 




[B(k)] 


b 

C 

C. 

1 

c 

D( ) 
[D] 

[DPHX 

[DX] 

[DZ] 

[EH] 

f 

g 

[H] 

[HAj] 

[HAW,] 

Hj(p)' 


[K]"^ r-i[M] + i p [A(lk)] 

wing span 
constant 

ratio between and (Section 5«2.3) 



reference chord 
flutter determinant 
viscous damping matrix 

matrix relating control system displacements to structural 
displacements 

differentiating matrix 

interpolation and differentiation matrix relating slopes 
at downwash collocation points to displacements at struc- 
tural nodes (equation 4.13) 

interpolation matrix relating translations at downwash 
collocation points to displacements at structural nodes 
(equation 4.13) 

extrapolating matrix 

frequency, Hz ‘ 

structural damping 

interpolation matrix relating displacements at lumped 
aerodynamic load points to displacements at structural 

nodes (equation 4. IT) 

j = 0,1, 2, 3; defines [HA(k)] = [H]'^[AIC(k)] in cubic spline 
interpolation (similar to equation 4.37) 

J = 0,1, 2, 3; defines [HAW(k] = [H]'^(AIC(k)l [W(k] in cubic 
spline interpolation (similar to equation 4.37) 

transfer function of automatic control system (function 
of p) 


VI 


h. 

1 

TD(x,y) 

[IPHX] 

[K] 

[Kq] 

[AK^] 

[K] 

k 

KEAS 

^^(k) 

[M] 


[M] 


[M,] 

[AM^] 

M 

m. 

1 

P. 

1 

P 


displacement at aerodynamic load point in deflection mode 

imaginary part of D(x,y) 
interpolating matrix 

stiffness matrix, modalized stiffness matrix 
base stiffness matrix 

incremental stiffness matrix per unit design variable i 
(l + ig)[K], stiffness matrix 

reduced frequency k = 
knots equivalent airspeeci 

polynomial multip].ier-5 used in Lagrange’s interpolation 
formula 

mass matrix, modalized mass matrix 
r 

LMJ , mass matrix 
c 

base mass matrix 

incremental mass matrix per unit design variable i. 

total mass associated with the design variables 

design variable associated with structural mass (in mass 
or weight units) 

general design variable 

complex root of flutter equation for a given flight 
condition (p = (Y+i)k) 

aerodynamic lifting pressure' distribution corresponding 
to deflecd.lou.fflbdd j ' 

vii - 




'A 


q 


modal degrees of freedom (modal participation coefficients) 

;i 

: 

{q} 


modal column corresponding to solution of characteristic 
equation 

■ SI 

ij 

'll 

RD(x,y) 


real part of D(x,y) . 

>! 

r 

H 


modal row corresponding to solution of characteristic 
equation 

i. 

t 

I 

T(x,y) 


defined by equation 2.20 

1 

i 

V 


speed, flutter speed 

I 

f 

II 

''r 


required flutter speed (minimum allowable flutter speed: 

j] 


1.20 Vp for commercial, 1.15 Vp for military) 

1 

vu 


unsatisfactory flutter speed 

1 

5 

''d 


design speed according to Federal Aviation Regulations 

g 

!i 

iii 

[VPIC] 


matrix of velocity potential influence coefficients 

i-f 

1 

"l 


total mass addition associated with the additive column 
of design variable increments (equation 5-21) 

n 

1 

1 — 

1 

s: 

ro 

ftJ- 

j^W^j row matrices for numerical integration 

■7 

1} 

ii 

1 

[w],[w(k>] 

= [DX] 

+ ik[DZ] angle-of-attack generating matrix (funcion of k) 

i 



(equation J^.l6) 

i 

[WF],[WFD] 


generalized weighting matrices (Section U.1,2) 

I 

X 


coordinate in a fore-and-aft direction, general variable 

1 

1 

y 


coordinate in a lateral direction, general variable 

H 

P 

1 

{Z} 


column matrix of aerodynamic forces 

1 

{z} 


column matrix of displacements of structural nodes 

1 

[t] 


matrix of modal columns of displacements of structural 
nodes - - - 

1 

\ 

a 


angle of attack 

j 

i 

i 



general design variable 

I 

] 


I 

I 

viii 




Y 


X 

X 

P 

<j> 


( 1 ) 

t 


T.E. 


normalized real part of p = (y+i) k 
(l+ig )¥'■', Lagrangian multiplier 
Lagrangian multipliers 

air density 

aerodynamic velocity potential 
circular frequency, rad/s 

indicates derivative of a one-variable function 
trailing edge 


in-flight mode: modal column corresponding to a characteristic solution of 

the flutter equation 


flutter mode: in-flight mode that becomes unstable within the velocity 

range considered 

hump mode: in-flight mode with a minimum damping point within the 

velocity range considered 

Identification of count formulas XYyZ (Section f+.3.2) 


X = Oper 
X = St or 
X = Cor 

X = Read 
Y 

y = a 

y = b 
Z = Poly 
Z - Spline 


total number of operations 

number of matrix element words to be stored for easy access 

number of matrix element words to be kept in core for 
interpolation 

number of matrix element words to be read into core for 
interpolation when switching from one k interval to the 
next 

sequence of multiplication options (equations 
[W(k^)] is input 

[DX] and [DZ] are input 
cubic polynomial interpolation 
cubic spline interpolation 


ix 


t 


I 

1 

I 

{ 



r^r" 

The interpretation of the matrix dimensions used is as follows: 

M: ntunher of modes used in flutter analysis 

N: number of discrete degrees of freedom: structural displacements 

K: number of lumped aerodynamic forces 

D: number of downwash collocation points 

<j>: fraction of non-zero elements in rows of [W(k)] and [H] 

. V ' 

In addition: 

n_: mimber of k intervals needed 

0 - 

nj + 3 




STUDY OP FLUTTER RELATED COMPUTATIONAL PROCEDURES FOR 
MINIMUM WEIGHT STRUCTURAL SIZING OF ADVANCED 
AIRCRAFT - SUPPLEMENTAL DATA 

By R. F. O’Connell, H. J. Hassig and N. A. Radovcich, 
Lockheed-California Company 


1. INTRODUCTION 


In a companion report (Reference l) the authors have presented the 
results of the present study that are in most direct response to the objec- 
tives of the contract. To generate the results in Reference 1, however, 
considerable supporting work was performed that has potential significance 
for the worker in the field of structural optimization. The present report 
summarizes a major portion of that supporting work. 

The two-dimensional Regula Falsi approach to solving the generalized 
flutter equation is presented in considerable detail including some experi- 
ence with a frequency and velocity sweep method for a global search for 
flutter frequency and speed. 

The method of Incremented Flutter Analysis and some of its applications 
is presented. 

In a section on aerodynamics a method is presented for the use of 
velocity potential aerodynamic influence coefficients in the basic five 
matrix product formulation of the generalized aerodynamic force coefficients. 

A method to prevent oscillatory non-convergence, "hunting", in the flutter 
solution due to cubic polynomial interpolation of the aerodynamics is pre- 
sented. Options for computational operations required to generate general- 
ized aerodynamic force coefficients are compared in detail. 

Section 5 is an attempt to approach flutter optimization with the help 
of elementary considerations. The authors believe this section can aid the 
worker in the field of flutter optimization in improving methods of analysis. 

Practical experience with an ad hoc method of optimization is presented. 
It indicates that the number Of modal degrees of freedom required for a m,od- 
erately complicated arrow wing supersonic transport design may be considerably 
higher than the momber used in the numerical examples found in the literature. 

• c: Expressions for the derivatives of , the in-flight modal damping with 
respect to the design variables were derived during the study. They are pre- 
sented, together with other derivatives, in Section 7. 

The report concludes with a bibliography. • - 
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2. TWO-DBffiNSIONAL REGULA FALS 


One of the objectives of the present study is to survey and evaluate 
methods of solving the flutter equation which are suitable for repetitive 
use in a structural resizing procedure with flutter constraints. Several 
methods have been evaluated and the results are presented in Reference 1. 

The two-dimensional Regula Falsi approach was chosen for further development 
since it was considered most promising after an overall engineering evalua- 
tion, including the generation of some numerical data. 

In the following sections the basic mathematical formulation of the 
method is developed and the generalized flutter equation to which it is 
applied is discussed. One section is devoted to freqiiency and velocity sweep, 
a global search method for flutter speed and frequency that was partly devel- 
oped as a result of two-dimensional Regula Falsi investigations. Finally the 
Lockheed-California Company's two-dimensional Regula Falsi program as devel- 
oped for use in a flutter optimization is discussed. 


2.1 Basic Mathematical Formulation 

The Regula Falsi and Newton's method are related iterative methods of 
solving the non-linear equation f(x) = 0. The basic Regula Falsi is shown 
in Figure 2-1. The function f(x) is evaluated at two trial values of x, 
xp and X 2 , leading to fp and fp. By linear interpolation or extrapo- 
lation a new value of x, x^, is found with the associated f 2 . Next f^ and 
f 2 are used to generate x^. In a variation of the method, here called 
Type II, the interpolation or extrapolation is done with the latest value of 
f and the smallest of all preceding f's (Figure 2-2) 

The basic Newton Method is shown in Figure 2-3. The function f(x) and 
its derivative f’(x) are evaluated at a trial value x^, leading to fp 
and fp .: The tangent to the curve at Xp is intersected with the abcissa 
leading to a new value of x, X 2 , after which the procedure is repeated. The 
derivative f'(x) can be determined analytically or by means of a finite 
difference technique. In the latter form Newton's methpd resembles the 
Regula Falsi method even more than in its basic form. 

Figure 2-4 illustrates how a poor initial trial xp may lead the iter- 
ation process away from the desired solution in the case of the basic Newton 
method. The finite difference’ form of Newton's method might have done better, 
especially if a relatively large value of Ax were initially chosen. 

Obviously the Regula Falsi also can have convergence troubles. However, 
ty trying to choose the values of xp and Xg such that they tend to 
straddle the solution these convergence troubles are minimized. 

Both the Newton method and the Regula Falsi can be expanded to n 
equations with n unknowns, but, even for n=2, a graphical representation 
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Find 


Find Xj fi*om f(xg) and 

Xj^ from fCx^) and f(xg), 
x^ from f(xj^) and f(x2), 
x^ from f(x^) and f(xj^). 


Figure 2-1: Regula Falsi - Type I. 


Find Xg from f(xg) and f(x^); 

xj^ from ^(x^) and smallest of f(xj^) and f(xg); 
x^ from f(xi^) and smalle.;* of f(x^),f(x2) and f(x2); 
xg from f(x^) and sm^est of f(xj^),f(x2),f(x2) and f(x^) 
(xg notli illustrated) 


Figure 2-2: Regula Falsi - Type II. 








of the proced-ures is almost impossible. During the present contract, it was 
found that for n = 2 the type of non-convergence suggested by Figure 2-k 
occurred rather frequently, even when the Regula Falsi was used. 

In the following the Regula Falsi is presented in an extension to two 
equations with two unknowns. In this form it is referred to as the two- 
dimensional Regula Falsi. . 

Let 


F(x,y) = 0 and G(x,y) = 0 


( 2 . 1 ) 


be the two equations with tile two unknowns x and y. Let (^"-^>2,3) 

be three pairs of trial values in the vicinity of the solution to the equa- 
tions 2.1. Define a plane f(x,y) = a^x + b^ + c^ containing the points 

F(x. , y. ) (i=l,2,3) and a plane g(x,y) =ax+by+c containing the 
XI g g g 

points G(x^,y^). 

The coefficients a^, b^ and c^ are determined by the matrix equation: 


^1 ^1 


X2 


^3 ^3 


r ^ 

a. 


y = < 




F(Xi,yi) 

F(x2,y2) 


F(x3>y3) 


{ 2 . 2 ) 


In an obvious shortened notation: 


r ^ 

a 

f 


‘'f 




[XY] 


r~l 


(2.3) 


Similarly, 


r > 

a 


r ^ 

G. 


-1 


b„\= [xyI^Vg 


g 




i2.h) 


5 


' Figure 2-5 illustrates the procedure. The solution to the 

equations f(x,y) = 0 and g(x,y) = 0 is taken as the first approximation 
to the equations 2.1. This leads to the equations 

VU 

a X, + h y, = - c (2 

g 4 g 4 g 

from which: 


l..l 


\ ‘’f' 

-1 

i 

•=f • 

P1.J 


a b 
_ g g . 


c 


(2.7) 


The process is repeated using the pairs xlj,,yi^ and two of the pairs 
Xi»yi (i=l,2,3), eliminating the pair that least satisfies equation 2.1 
(e.g. the pair with the largest values of |F(x^,y^)| + |G(x. ,y^)|). It 

is continued until convergence is reached. 

Characteristic of the two-dimensional Regula Falsi is: 

1. three pairs of trial values are required to initiate the process 

2. the two functions, but no derivatives, must be evaluated once at 
each step, except at the first step for which the two functions 
must be evaluated for three pairs of values x,y 


2.2 The Generalized Flutter Equation 
When using the k method the fluter equation can be written as: 

i [M] k^ + [K] - I p [A(ik)]j {q} = 0 (2.8) 

One of several possible methods of solving this equation is to deter- 

1+ig 

mine the characteristic value X = - a -— for several values of the reduced 

V 

frequency k = ^ , keeping all other quantities in the equation constant 
(Reference 2). 
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In the p-k method the flutter equation is 


[M] + (1+ig) [K] - I pV^ [A(ik)j 


{q} = 0 


(2.9) 


and solutions p=(Y+i)k are sought for selected combinations of values of 
V and p, such that [A(ik] is evaluated for the value of k that defines 
the imaginary part of the solution. The p-k method formulation is con- 
venient for the inclusion of viscous damping and control system transfer 
functiohti (Reference 3). This is accomplished by writing: 


-^[M]p^ + ^D]p + (l+ig)[K] -|pV^[A(ik) - IH.(p)[Dj] 
c 


{q} = 0 


( 2 . 10 ) 


where H.(p), J=l,2.., represents transfer functions of the control system 

^ • 

and [D.] relates control system displacements to the structural displace- 
ments; [D] is a viscous daittiping matrix. 


A further generalization of the flutter equation can be made by making 
the stiffness matrix and the inertia matrix functions of design variables 

which is the standard procedure for structural optimization. In addi- 


^i » 

tion, other quantities, such as [D] , [D.J, as well as transfer function 

J 

coefficients in H.(p) may be made functions of design variables. 

Equation 2.10 implies that the determinant of the square matrix on the 
left hand side ((is zero and thus, in a very general form, the flutter equation 
can be written Vas : 




D {(Y+i)k,g,V,p,m } = 0 


( 2 . 11 ) 


The quantity D is called the flutter determinant. For arbitrary valuep of 
the variables it has a complex Value. Thus equation 2. Xl represents ttiro 
equations equivalent to equations 2.1. In principle, they can be solved for 
two unknowns for given values of the other variables. 


Letting y ~ 0 solving for g and V corresponds to the tradi- 

tional k method of solving the flutter equation. Equation 2.8 can be 


written in the form 


[ 


X[lJ - [B(k)] {q} = 0 


] 


( 2 . 12 ) 
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by the classical power method and improve- 


1+is 

and can be solved for X =' — — 
ments thereof. 

Solving for (y+i)k corresponds to the p-k method. According to 
Reference 3, equation 2.10 for (small) values of y of interest is "almost** 
an analytic function of p = (y+i)k and it can be solved by determinant 
iteration which is basically a one-dimensional Regula Falsi method. 

In the k method and the p-k method the flutter speed is determined 
indirectly by interpolation between solutions for several values of k or 
several values of V. However, equation 2.11 can be solved directly for the 
flutter speed at a given value of the structural damping g by letting 

Y = 0 and solving for the two unknowns k and V. 

Solving equation 2.11 for k and one of the design variables, assuming 
all other variables fixed, has found application in structural optimization 
with flutter constraints (Reference l). It forms the basis for the method 
of Incremented Flutter Analysis which is discussed in Section 3. 

Since equation 2.10 actually is not an analytic function of p - (y+i)k 
it should be considered as two equations with the two independent unknowns 

Y and k. Lockheed's program for solving the flutter equation according to 
the p-k method has indeed an option to use the two-dimensional Regula Falsi. 
During numerical testing it was found to be slightly less efficient than 
determinant iteration in that on the average one more iteration step is 
required. 

The two-dimensional Regula Falsi Program used in conjunction with the 
p-k method is a general purpose multi-dimensional Regula Falsi program. 
Initially it was used unchanged to solve the generalized flutter equation for 
k and V, and for k and . Few convergence difficulties were encountered 
during an optimization study of an arrow wing supersonic transport when solv- 
ing for k and mj_. Initially, however, convergence difficulties occurred 
erratically when solving directly for flutter speed and flutter frequency. 

To gain an understanding of the problem the determinant al flutter equation 
for several fixed configurations was written as 

D(a),V) = 0 ■ (2.13) 

and its behavior as a function of m and V was studied in depth. This 
led to an improved version of the two-dimensional Regula Falsi program, 
described in Section 2.U, and had as a by-product a potential globsl search 
method for estimating flutter speed and flutter frequency which is discussed 
in the next section. 
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2.3 Frequency and Velocity Sweep 

As indicated in the previous section initial numerical experiments with 
the two-dimensional Regula Falsi ran into convergence difficulties when solv- 
ing directly for flutter speed and frequency. To gain an, understanding of 
the problem the flutter determinant, D((jo,V), for a fixed configuration was 
evaluated for many values of o) and V. Various types of plots relating w, 

V and the real and imaginary part of the determinant were studied. One type 
of plot, the argument (or phase) of the flutter determinant versus frequency 
and versus velocity, led to results that may be of value in developing a global 
search for approximate values for the flutter speed and flutter frequency. 
These approximate values could subsequently be used in a local search by the 
-two-dimensional Regula Falsi method to determine accurate solutions . Since 
in most approaches to flutter optimization it is assumed that approximate 
solutions of the flutter equation are available to initiate a new solution 
process this global search method was not pursued in the present study. The 
results obtained thus far are summarized in this section. 

2.3.1 Frequency Sweep . It was found that if the complex flutter determinant 
is evaluated for a speed V^^ which is less than the lowest flutter speed the 

values of « that satisfy the equation arg {D(u,V^)} = ^ are good approxi- 
mations for the flutter frequencies. Three sample cases are presented and 
discussed. 

All three cases are for symmetric flutter analyses of the arrow wing 
transport shown in Figure 6-1 at a Mach number M = 0.9 and an operating 
weight empty (OWE) of 321,000 lbs. Further details regarding structural rep- 
resentation, structural and vibrational degrees of freedom, and aerodynamics 
are presented in Section 6.1. For all cases the structural damping g — 0.02. 
Natural vibration modes, including two zero frequency modes, are used in the 
flutter analysis. 

Figure 2-6 is the f-g-V diagram for Case 1. It results from an 
analysis with 20 modes and variable air density (matched atmospheric 
conditions). The flutter frequencies and associated speeds are 2.07 Hz, 

318 KEAS and 2.hh Hz, 375 KEAS. Figure 2-7 shows the argument of the flutter 
determinant as a function of the frequency for several values of V. It can 
be seen that only those flutter frequencies are detected that correspond to 
a flutter speed above the value of V for which the argument is determined. 

Figure 2—8 is the f-g— V diagram for Case 2. It results from an 
analysis with 20 modes and a constant air density corresponding to 18,000 ft. 
To confirm the true character of the fourth and the fifth mode, solutions for 
these modes between 400 and 420 KEAS were determined at very small velocity 
intervals (Figure 2-9). The flutter frequencies and associated flutter speeds 
are 1.15 Hz, 542 KEAS; 2.71 Hz, 390 KEAS and 4.21 Hz^ 609 KEAS. Figure 2-10 
shows the frequency sweep: argument of flutter determinant versus frequency 
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Figure 2-7: '^ase 1; Frequency Sweep for Several Speeds; f is Flutter 

Frequency Obtained from Figure 2-U. 
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'-10: Case 2; Frequency Sweep for 250 KEAS; is Flutter 

Frequency Obtained from Figure 2-8. 





for V = 250 KEAS. The three flutter frequencies are detected. Another 
flutter frequency f » I .70 Hz is indicated. It is believed if the f-g-V 
diagram were extended to higher velocities a flutter speed corresponding to 
1.70 Hz would be found. 

The f-g,-V diagram for a third case is shown in Figure 2-11. This case 
is different, from Case 2 in that only 12 modes are used. This f-g-V diagram 
is very similar to the one in Figure 2-8. The key difference is that the 
curve corresponding to the lowest flutter speed in Figure 2-8 has an S-curve 
above the zero damping line and continues steeply into the unstable region, 
whereas in Figure 2-11 the corresponding curve turns stable and becomes a 
hump mode. 

The zero damping points in Figure 2-11 occur at 2 . 7 I Hz and 4l0 KEAS, 

3.05 Hz and 4^0 KEAS, l.lU Hz and 552 KEAS and, on the back side of the hump 
mode, 2.85 Hz and 451 KEAS. 

The frequency sweep corresponding to Figure 2-11 is shown in Figure; 2-12; 
it is very similar to Figure 2-10. The flutter frequencies corresponding to 
the backside of the hump mode and the steeply unstable mode, 3-05 Hz and 
2.85 Hz respectively, are not detected. Additional cases very similar to 
Case 3 were investigated with similar results. 

A preliminary conclusion is that the f 4quency sweep as defined above 
can be used for a global search for flutter frequencies. However, there is 
no assurance that all flutter frequencies are detected. Further investiga- 
tion may lead to a useful method. 

2 . 3.2 Velocity Sweep . It was found that if the complex flutter determinant 
is evaluated for a frequency ,o)j_ which is a 'little less than a flutter fre- 
quency the value V that satisfies the equation arg {D(w^,V)} ~ ^ 

is a good approximation for the associated flutter speed. 

Figure 2-13 shows the argument of the flutter determinant for Case 1 as 
a function of the velocity for several values of f . The highest /alue ,.of f 
that is still below the estimated flutter frequency according to Figure 2-7 
is f = 2.00 Hz. For the corresponding value of w, arg {D(w, ,V)} = - is 

satisfied for V = 319 KEAS. The actual fluttex* speed is 318 KEAS (Fig- 
ure 2-6 ) ., 

A velocity sweep for Case 2 is shown in Figure 2-l4 for f = 2.65 Hz, 
a little below the estimated flutter frequency of 2.73 Hz indicated in 
Figure 2-10. No intersection with the -90*^ argument line occurs. A velocity 
sweep for Case 3, for f = 2.70 Hz, is shown in Figure 2-15; f = 2.70 Hz is 
a little below the estimated flutter frequency of ,2.77 Hz indicated in Fig- 
ure 2-12. The intersection with the -90® argument line indicates a flutter 
speed of 4l2 KEAS which compares with 4l0 KEAS as obtained from the f-g-V 
diagram (Figure 2-11). 
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Figtire 2-12: Case 3; Frequency Sweep for 250 KEAS; is Flutter 

Frequency Obtained from Figure 2-11. 
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Figure 2-l4: Case 2; Velocity Sweep for f = 2.65 Hz; V_ is Flutter 

Velocity Obtained from Figure 2-8. 
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It is concluded that the velocity sweep is an uncertain method for a 
global search for flutter speeds. Possibly further investigations may lead 
to an effective approach. 

2 . 3.3 Conclusion . The frequency sweep and velocity sweep, if made reliable, 
would provide a global search method for flutter speed and frequency that 
does not require the solution of characteristic value problems and that is 
free from the problems of mode following associated with the usual modal solu- 
tions. This potential gain justifies further examination of this approach. 


2.U Description of the Basic Two-Dimensional 
Regula Falsi Program 

Lockheed-California Company's two-dimensional Regula Falsi program 
solves the equation D(x,y) = 0 for x and y, given a reasonable initial 
estimate of x and y. 

The two equations implied by D(x,y) = 0 are written as: 

RE>(x,y) = 0 and ID(x,y) = 0 (2.15) 

The two-dimensional Regula Falsi i^iprocedure is a separate sub-program 
which requires the values of RD(x,y) \,and ID(x,y) to be supplied through 
a calling list. RD(x,y) and ID(x,y)\are therefore generated outside the 
Regula Falsi sub-program and can easily%e adjusted to the problem at hand. 

The Iteration process is initiated/by three pairs of values of x and 
y as described in Section 2.1. These ^Jyalues are defined as follows: ,, 

i' 

^1 " ^0 ^ ^1 "" ^0 

^2 = ^1 ^I ’ ^2 = ^1 ^ 2 . 17 ) 

X 3 = Xi ; y 3 = yi + R^ Ay ( 2 , 18 ) 

Xq and y-Q are estimated values of the solution, defined by the User. For 
the definition of the second and third trial pair the user must define Rt, 

Ax and Ay. 

Two modes of iteration are provided: Mode US (unrestricted stepsize) 

and Mode RS (restricted stepsize), to be defined later. A complete iteration 
may consist of one or more searches in Mode RS onljr, one search in Mode US, 
or combinations of Modes US and RS. If a combination is used. Mode US 
searches are conducted first, followed by the same number, or one less, of 
Mods RS searches. The user specifies the type of search desired and he 


22 


controls the number of searches by defining the maxlmiun number of iterations 
per search and the maximum total number of iterations. Each time a new 
search is entered the iteration is reinitiated by defining new values of 
Xi,yi, equal to the values of the search just terminated. The 

quantities Xjjj and yj^^ are defined in connection with equation 2.21. 


The quantities Ax 
iteration the values of 


and Ay define a search region. 
X and y are restricted to: 


During the entire 


Xq - Ax < X < Xq + Ax 
yQ “ Ay < y < Yq + Ay 



(2.19) 


This restriction is necessary to prevent the iteration p,rocess from 
being led away from the desired solution. If during the iteration either x, 
or y, or both exceed the boundary, the corresponding x or y that is in 
violation of the boundary is replaced according to a formula that depends on 
the mode of iteration, as discussed in a subsequent paragraph. Note that the 
boundaries are based on Xq and yQ and do not change when a new search is 
initiated. 


The Type II Regula Falsi is used (Figure 2-2). Assiime that the trials 
x. ,y., y. have led to the next iterated value x,^^, y,. 


i+l» •^1+1* i+2* •'i+2 

Then for the next iteration step ^1+3 


i+3* •^i+3* 
is used together with those values 


of x.,y.(j<i+3) that correspond to the smallest two values of the quantity 
J J 

T(x,,y.) = T„(x,,y,) + TT.(x,,y,) which is defined by equation 2.20. 

J O ^JJ -i-JJ 


T (x.,y ) = the larger of 
b 0 J 


LD(Xj,y^) 


LD(x^,y^) 


-lO”^ and 0; L = R,I 


( 2 . 20 ) 


where lO”^ is a convergence criterion (equation 2.25). 


For this purpose the program retains the values of RD and ID correspond- 
ing to the smallest two values of T(x^,y^), and the associated »yj pairs as 

they occur during one search. Whenever a smaller T is encountered it replaces 
the larger of the two that were retained. When a new search is entered the 
x,y corresponding to the smallest T of the previous search become the x^ jy^ 

of the new search. 


The distinction between the two modes of iteration, RS and US, lies in 
the stepsize. Stepsize is defined by 



X. . o 

1+3 


- x 


m 


6y = I yi+3 - 


( 2 . 21 ) 


where x correspond to the smallest value of T during the current search, 

m’ m 
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In Mode US (unrestricted stepsize) there is no restriction on 6x and 
6y except that eq.uation 2.19 must he satisfied. In mode RS (restricted 
stepsize) 6x and 6y are restricted according to the equations: 

6x < 5y < (2 


The user defines the value of Rg. 

When equations 2.19 and/or equations 2.22 are violated, corrective 
measures are taken. Several possibilities are recognized and treated as 
follows : 

1. Mode RS and one of the equations 2.22 is violated; e.g. 

6x > R R^Ax. Then x resulting from the last Regula Falsi 
step is replaced by 


x.^_ = X + 

1+3 m 


i+3 


m 


l^i+3 " 


R R.j.Ax 

D X 


( 2 . 23 ) 


Similarly for y. 

2. Mode RS and one of the equations 2.19 is violated. The quantity 

\ 

in violation, say is replaced by ^^+3 according to equa- 

tion 2 . 23 . If x ^+3 still violates equations 2.19, the iteration 
is terminated. 

3 . Mode US and one of equations 19 is violated. The quantity in 

violation, say i^ replaced by x ^^3 



m 


+ 


./i±3 


X 


i+3 



2R^Ax 


(2.24) 


If ^j_ 4_3 still violates equations 2.19, the iteration is 
terminated. 
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There are three further conditions under which the iteration is 
terminated: 

1. A converged soi.ution is obtained. The convergence criterion is: 


^converged ^ ^ 

“converged < 


(2.25) 


where e -is defined by the user. 

2. There are two consecutive occurrences, or three occurrences in total, 
of either x or y exceeding the boundaries defined by equa- ■ 

tlon 2.19. 

3. The number of iterations exceeds the total number of Iterations 
permitted as specified by the user. 

If the iteration is terminated without convergence either the entire 
calculation is terminated, or a recovery module is entered. At present there 
is one recovery module, called the Z-module, which is described in the next 
section. 


2 . 5 The Z-Module 

The Z-module provides a recovery procedure for the two-dimensional 
Regula Falsi program if the iteration is terminated without convergence having 
occurred. It is specifically designed to provide recovery under the circum- 
stances encountered in an automatic resizing program with flutter constraints. 
It can be used, however, where the solution for a basic configuration is known 
and the solution for a modified solution is sought. 

Let X stand for any combination of the variables x and y. Let 


D (X) = 0 be the determinantal equation for the base configuration, and 
B 

Dg(X) + AD(X) = 0 the equation for the modified configuration; X^ is the 

known solution of D_(X) = 0 and X is the estimated solution of 

B U 

D^(X) + AD(X) =0. 

Assume that solving Dg(X) + AD(X) =0 with X^ initiating the itera- 
tion led to a termination without convergence. The Z-module then generates 
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a new determinantal equation DB(X) + Z AD(X) = 0 and first trial value 
X = X_ + Z(X_ - X^). Z is a fraction smaller than unity and in effect the 
Z-module reduces the difference> between the modified and the base configura- 
tion. If X„ is chosen equal to X , a convenient choice in flutter optimiza- 

tion, X = X tends to be a better first trial value for D_(X) + ZAD(X) = 0 
0 B -D 

than for D^(X) + AD(X) = 0. 

The first value of Z is Z = 0.5. If this does not lead to convergence 

Z = 0.25 is used. In general Z = (0.5)^, where n is the number of succesr 

sive occurrences of non-convergence. The value of Z at which convergence 

occurs is identified as Z and the corresponding solution as X„. New trial 

values for the original problem D (X) + AD(X) = 0 are obtained by linear 

B 

extrapolation using X^, and X^. 

B U 

If non- convergence occurs again, the design condition at Zq becomes 
the new base design: 

Dg(X) = Db(X) + Zq AD(X) (2.27) 

and ^(X) = (l-Z ) AD (2.28) 

and n is set at n = 1 and Z is reset to 0.5., 

The reduction procedure is repeated until another value of Z at which 
convergence occurs is found. At this point the three known solutions are 
used in obtaining new triaT values for Z = 1.0. 

The above iterative process is continued until either convergence is 

obtained at Z = 1.0 or until n = wh,ere 

^MAX ~ allowable successive convergence failures. 


2.6 Numerical Examples 

Many niunerical examples have been used to develop and test the two- 
dimensional Regula Falsi Program and the Z-module. The example in table 2-1 
is chosen to illustrate the program. The final solution can be seen to lie 
outside the search region. The initial trial" is' made equal to the base 
solution. The first search is terminated because of excessive violations of 
the search limits. The Z module sets Z = 0.5, cutting the difference between 
the base problem and the new problem in half. When again thbre are more than 
the allowed number of search limit violations Z is reduced ho 'G ;-25 tfor 
which value convergence occurs. This solution and the base solution are^ used 
to generate an improved initial trial value for the new problem (Z=l) which 
leads to a converged solution. 


3. INCREMENTED FLUTTER ANALYSIS 

Incremented Flutter Analysis is a useful tool that can enhance several 
methods of optimization with flutter constraints (Reference l). It is a 
method for directly solving for the value , of a design variable that satisfies 
a given flutter constraint. In the original presentation of the method 
(Reference k) this basic idea was coupled with a formulation that led to a 
large reduction of the order of the characteristic value problem to be solved 
without losing any of the accuracy implied by the original order. This 
formulation is useful for problems with a, small number of design variables, 
such as the determination of external store parameters satisfying the flutter 
speed requirements. When the number of design variables is large, and when 
they can not easily be isolated in the mathematical representation of the 
problem (e.g., if problem is modalized) the reduction of the order of the 
characteristic problem to be solved becomes impracticable. Such is often the 
case with flutter optimization. In this section the formulation of 
IncrerAented Flutter Analysis for general application in flutter optimization 
programs is presented, together with applications. 

It should be noted that the essence of Incremented Flutter Analysis was 
recognized elsewhere (Reference 5). 


3-1 Form of the Flutter Equation 

The characteristic flutter equation can be written as 

D(g,Y,k,V,ra^) = 0 (3.1) 
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INPUT: Allowab3.e Iterations per search = 7, Allowable Iterations per Z = 10 

Allowable Iterations per run = 40, Xq = 2.89, Yq = ^53, Ax = .2, Ay = 20 

= .2, Rg = .5, e = lo"^, Xg = 2.89, yg = i+53 


RESULTS: 


Iteration 

X 

y 

RD 

— 

ID 

z 

r~ 

COMMENTS 

Number 






1 

2.890 

45 3. 00 

-0.8585 

i.i415 

1.00 


Initial Trials: ^1 ” ^ 0 ’ ^1 ~ ^0 

2 

2.930 

453.00 

-1.1542 

1.3529 



Xj, = Xq . R^AX, 

, 3 

2.890 

457.00 

-0.8343 

1.0854 



^3 “ ’'o’ ^3 ■ ^0 ^ 

k 

2.818 

465.00 

- 0.5746 

0.5006 



Upper y search limit violated. 

5 

2.738 

457.00 

-0.4335 

0.1322 



Lower x and y search-.J-imits violated. 
(Discontinued due to excessive search 
limit violations) 

5 

2.890 

453,00 

-0.8852 

1.1148 

0.50 


Initial trials are the same as at 

7 

2.930 

453.00 

-1.0788 

1.4557 



■ Z = 1.0 since x^ = Xg and y^ = Vg 

8 ■ 

2.890 

457.00 

-0.8186 

0.9703 



9 

2 . 810 

449.00 

-0.7219 

o. 4 i 4 i 



Lower x and y search limits 
violated 

10 

3.076 

457.00 

-1.9066 

0.8490 



Upper y search limit violated. 
(Discontinued due to excessive search 








limit violations ) 


Table 2-1; Numerical Example of Two-Dimensional Regula Falsi Iteration Including 
the Use of the Z Module 



RESULTS 


Iteration.: 

Number 

X 

11 

2.890 

12 

2.930 

13 

2.890 

l 4 

2,926 

15 

2.906 

16 

2.904 

17 

' 2.905 


U 53 .QO 

453.00 

457.00 

465.00 
4 t 0 . 6 T 
469 . 09 
469.32 


51a. 27 

518.27 

522.27 
517.68 

517.84 

517.88 


-0.8511 

-0.8566 

-0.6953 

-0.0107 

0.0958 

-0.0197 

0.0007 


-1.2027 

1.5326 

-0.5506 

-0.0542 

■0.0055 

0.00008 


0.1110 

0.0128 

0.007 


■0.7973 

0.4762 

-2.0998 

■0.0586 

0.0171 

0.0004 


COMMENTS 

Initial trials 

are the 

same 

as at 

Z = 1.0 since 

^0 = 

and 

II 

0 

Upper y search limit 

violated. 

Converged for Z 

; = 0.25. 



] X and y obtained by linear 

extra- 

1 polation using 


and : 

^7’ ^17 

Converged at Z 

= 1.00 




Table 2-1; 


Numerical Example of Two-Dimensional Regula Falsi Iteration Including 
the Use of the Z Module (Continued) 























The quantity D is called the flutter deterininant and has a real and 
imaginary part. Thus equation 3.1 represents' two equations. For this 
discussion D is considered a function of the structural damping g, the 
decay rate y > the reduced frequency k, the velocity V and the design 
var i able s m^ ( i =l^n ) . 

Since equation 3.1 represents two equations it can be solved for two 

unknowns. In the traditional k method of solving the flutter equation 

Y = 0, m. is fixed and for a series of values of k equation 3.1 is solved 
1 

for g and V. In the p-k method of solving the flutter equation g and m;j_ 
are fixed, and for a series of values of V the equation is solved for y 
and k. 


Letting y = 0 and giving g and mj_ fixed values, equation 3*1 can 
be used to solve directly for the flutter speed and the associated reduced 
frequency. This use of the flutter equation has gained importance in 
connection with structural optimization with flutter constraints. 

Incremented Flutter Analysis is characterized by the addition of m^^ 
to the variables of which D is a function. Thus the meaning and use 
of the flutter equation is generalized. When in equation 3.1 the values of 
g, Y, V and all but one of the m^^'s are fixed it can be solved for the 
value of k and the variable m^j^. If a real solution exists, the value of the 
unknown mi is found that, together with the fixed mj^'s defines a structiire 
that at the given speed and structural damping exhibits the given rate of 
decay. If y = 0, the given V is the flutter speed, and thus a structure 
has been defined with a preset flutter speed. 


3.2 Application to Finite-Difference Technique 


Most methods of structural optimization with flutter constraints require 
the evaluation of the derivatives of the flutter properties with respect to 
the design variables. Analytic expressions for the flutter-speed derivative 
(Reference 6) and the derivative of the logarithmic increment (Section 7*2) 
have been derived and applied successfully. However, derivatives obtained by 
the finite-difference technique have also been used successfully 
(Reference 7). This section shows how the forward finite difference quotient 
of the flutter speed with respect to a design variable can be determined with 
the help of the Incremented Flutter Analysis Technique. 


Let I m^ I define a structure with a known flutter speed V. Define a 

perturbed structure by » Where is obtained from { | i’y 

increasing the single j element by Am.. Assume this change increases the 

J 

flutter speed ah amount AV. 
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Write the flutter equation as: ' 

Y» k, V, {m^}, Am^ j = 0 (3.2) 

Let y - 0, g - structural damping, V the original flutter speed and 

{ } constant. Solve for each (i = l-^-n) separately, and the associated 

k value (which is of interest only for checking purposes). Each Am^ found 

is the amount that must be added to the i 'design variable in the perturbed 
system such that each Am^ by itself restores the flutter speed to its 

original value. By definition 


(Am.). . = - Am. 

1 i=J J 


(3.3) 


The values of Am^ thus determined can be used to approximate derivatives 
of the flutter speed with respect to the design variablejS. 

aV AV , , , 

am, ~ Am, ^ 

l x 

This requires the numerical evaluation of AV. For certain procedures, only 
normalized flutter speed derivatives are needed. Let the subscript R refer 
to a reference design variable. Then: 


av 

3m^ AiOp 

av Am, 

'“e 

and there is no need to compute AV. 


3.3 Determining the Magnitude of a Resizing Column 

When raising the flutter speed from an initial, deficient value to the 
required value, and in several methods of structural optimization with 
flutter constraints, it is necessary to determine the magnitude of a given 
design variable distribution column such that the flutter equation is 
satisfied at a given flutter speed. 

Incremented Flutter Analysis provides a technique for determining this 
magnitfUde directly, without need of generating several solutions to be used 
in an interpolation procedure. 



norm 
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The statement of the problem is : given a structure with an unknown 

flutter speed, defined by » and a distribution of design variable 

increments {Am^}, determine a scalar C such that the structure defined by 
{™i } "*" ^ ^ given flutter speed V. 

The solution is found by writing equation 3.1 as: 

D (g, Y, k, V, C|Am^}) = 0 ( 3 . 6 ) 

and solving for C and k. 

An extension of this procedure can be used to determine the relative 
magnitude of two resizing col\amns as might be needed in a one-dimensional 
minimization routine such as the one described in Section 6.6 of Reference 1. 
For this purpose equation 3.6 is written in the form: 

D(.g, Y> k, V, C C {Am^}) = 0 (3-7) 

Let Y = 0, and g, V, J m^ J , { Am. } and { Am. ) be fixed. For 

# 

successive values of C the corresponding values of C can be determined 
such that the flutter speed has the desired value. The one-dimensional mini- 
mi z at ijDn_t hen consists of determining the value of C .for which 
AM = CEAm. + CZAm. is minimum. 

11, 

3 . h Concluding Remarks 

The finite difference technique using Incremented Flutter Analysis has 
has been used successfully for approximating flutter speed derivatives in 
the numerical evaluation performed for the method of optimization described 
in Section 6.6 of Reference 1. It was concluded, however, that as the nimiber 
of design variables with respect to which the derivatives need to be 
determined increases, the analytical approach to determining flutter speed 
derivatives (Reference 6) becomes relatively more efficient, nevertheless, 
the technique described may be useful in special cases. 

Determining the magnitude of a column of design variable increments such 
that a flutter speed is exactly satisfied' can be a useful tool in all methods 
of optimization discus^:ed in Reference 1, except the penalty function method. 
It eliminates the need for more than one solution of the characteristic 
flutter equation and subsequent interpolation. 
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4. AERODYNAMICS 


One of the objectives of this study is to determine general, efficient 
and accurate computational forms in which to cast the unsteady aerodynamic 
parameters necessary for use in structural resizing required to meet flutter 
constraints. This effort, however, did not involve evaluation of aerodynamic 
theories. A general discussion and conclusions are presented in Reference 1. 
Additional discussion is presented in this Section. 

Reference 1 indicates that an aerodynamic theory based on the use of the 
aerodynamic velocity potential leads to the same basic five matrix product 
for the generalized aerodynamic force coefficients as the kernel function 
approach of Reference 8. Section 4.1 presents more detail. 

The problem of "hunting" when using aerodynamic matrices, interpolated 

for arbitrary values of k = may occur if the interpolation is piecewise, 

i.e., each interpolation covers only part of the range of k values. This 
problem is mentioned in Reference 1 as related to jumps in the derivative 
when the third degree polynomial interpolation formula is simply differ- 
siibiated. Section 4.2 discusses a form of "hunting" which can occur even if 
no derivatives are used. 

As part of the present study an effort was made to come to definitive 
conclusions regarding details of an efficient procedure for computing the 
generalized aerodynamic force coefficients for use in the repetitive calcu- 
lations required in a structural weight minimization procedure. Formulas 
for the number of computational operations needed to generate the aero- 
dynamics matrix have been derived and sample calculations have been made. 

This work is summarized in Section 4.3. 


4.1 The Use of Velocity Potential Influence Coefficients 

A generally useful formulation of the unsteady aerodynamic forces is one 
in which a matrix of aerodynamic influence coefficients,, [AlC(k)], expresses 
liamped aerodynamic loads at an aerodynamic loads grid in terms of angles of 
attack at downwash collocation points: 

{Z} = [AIG(k)] {a} (4.1 

Such a matrix can be obtained regardless of which aerodynamic theory is used. 

Two types of theoretical approach are widely used: one is based on the 

use of the acceleration potential, the other on the use of the velocity 
potential. Three methods of computing unsteady aerodynamic forces are 



conimonly used: subsonically , the integral kernel formulation of Reference 8 

and the finite element approach of Reference 9, both based on the use of 
the acceleration potential; supersonically, the supersonic mach box approach 
of Reference 10, which is based on the use of the velocity potential. 


Recently the velocity potential approach has received increased attention 
as a common basis for the computation of subsonic, transonic and supersonic 
aerodynamics (References 11 through 15), and the question arose whether 
a formulation in terms of velocity potential influence coefficients, VPIC, 
might have computations’! advantages when used in an optimization procedure. 
Therefore in the following, theoretical expressions are developed to express 
the generalized aerodynamic force coefficients in terms of the five matrix 
product of Reference 1, assinning velocity potential Influence coefficients 
are available. 

4.1.1 Basic equations. The generalized aerodynamic force is defined by 




h^dxdy 


(H.2) 


where pj is the lifting pressiire distribution associated with the j mode 
and h^ the displacement as socia^bW with the i"^^ mode. 


The pressure distribution can be expressed in terms of the velocity 
potential <j). The linearized Bernoulli equation gives: 


= -2 (If - m) 


(1|.3) 


where k is the reduced frequency: k = — • 

Thus equation h.2 becomes: 


A. . 

ij 


-4 




) “i 


dxdy 


(J+.4) 


Direct integration of equation' U.4 gives: 

/•T.E. 


A. . = • -2 


/•b/2 

f 

H 


J-b/2 



ik(|.j l^^dx + 


L.E. 


J ^T.E. 

1 9x 

L.E. 


dx 




(4.5) 


Integrating the term with by parts gives: 


A. . = -2 

ij 


/•b/2 ^T.E. 

I IHj h.dx ^ (h.*j 

J-b/2 Jl.S. 

J rT.E. 

L.E. 


E. 


(it. 6) 


8h. 

♦j sr 


In performing the integrations in equations k.^ and 1+.6 by niimerical 
integration one may consider first doing the chordwise integration at 
several stations y = constant along the span, and then doing the spanwise 
integration. In following this procedure it is found that inconvenient matrix 
formulations result in which mode-dependent and mode-independent parts are 
intermixed. Therefore the numerical integration is formulated as an area 
intregation. Thus equation I +.5 is written as: 


/•b/2 ^T.E. ^b/2 ^T.E. 

I h^dxdy - 2 I I h^ 

2 Jl.E. J-b/2 Jl.E. 


A. . = -2 

iJ 


J-b/2 

and equation k.6 as: 

•b/2 /»T.E. 


3d). 

J 

3x 


dxdy (4.7) 


^b/2 ^T.E. ^b/2 ^T.E, 

Aij = -2 I I ilc(j)^ h^dxdy + 2 j j 

J-bi2 Jl.E. J-b/2 Jl 


3h. 

♦j JT 


.E. 


-2 


/•b/2 

' 

J-b/2 


(4.8) 




Equation 4.8 has one more term than equation 4.7 but has the advantage 
that the integrand does not become infinite at the leading edge. 


^•1-2 The two-term integration . Written in a form for numerical integration 
equation 4.7 becomes: 


A. , = -2ik 
ij 


NWi 


•li - “ 


WHPI 


(4.9) 
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where L«ij and L«3j are integrating matrices associated with integration 

points at which h^ is defined. The option that j^^lj ^ included 

to provide the possibility of adjusting the integrating matrices to the shape 
of the functions to be integrated. 


Let hj^ be defined at y=constant stations at which (j)^ is defined, 
but assume the general case that hj^ and 6^ are defined at different 
points in the x direction. 


Let [IPHX] 
matrix such <|> . 
respectively, by 


be an interpolating matrix and [DPHX] a differentiation 
3(j). 

and — — ^ at the integration points are defined. 


[IPHX] {<!).} and [DPHX] {(f).} (^.10) 

J J 

Substituting the expressions i+.lO in equation 4.9 and interchanging the 
W and h matrices gives: 


A, . 

ij 



[IPHX] -2 



[DPHX] 




(ii.ii) 


which, with 
to equation 



[IPHX] = [WF] 


5.11 of Reference 1, 


and 


[«2] 


[DPHX] = [WFD] 


generalized to non-coinciding 


, corresponds 

h . and (b . 

1 1 


stations. The velocity potential can be expressed as the product of a matrix 
of velocity potential influence coefficients and an angle-of-attack 
distribution. 


{(f).} = [VPIC] {a } (4.12) 

J J 

The angle-of-attack distribution can be related to the structural displace- 
ments, {z}, by a differentiating matrix [DX] , generating , and an 
interpolating matrix ik[DZ], generating the frequency dependent equivalent 




[ 


[PX] + ik [DZ] 






(i*a3) 
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Combining the equations ii.ll, k .12 and U.13: 




IJ 


-2ik, J [IPHX] - 2 j^wj [DPHX] [VPIC] j^[DX] + ik [DZ]J {z^} 


Let [AIC] = 


-21k [wj 




[IPHX] - 2 [DPHX] 


[W] = [DX] + ik [DZ] 


[VPIC] (1+.15) 

{h.l6) 


[hj ‘ 


. I is obtained from the structural displacement z^ by interpolation: 


^h.J = l^z.J [nf (It.lT) 

With equations 4.15, 4.l6 and 4.17 equation 4.l4 becomes: 

A. . = I z. I [nf [AIC] [W] {z.} (4.18) 

ij L J 

For several modes z. and z. this becomes 

X j 

[A ] = [z]'^ [H]'^ [AIC] [W] m (4.19) 

1 J 

This equation is identical to equation 5*12 in Reference 1 and the [AIC] 
matrix can be evaluated by any suitable method. 

4.1.3 The three-term integration. Write equation 4. 8, .in a form for numerical 
integration: 


A. . - -2ik 
ij 


L’'iJ M 


} + 2 


W, 


8h. 
-1 


L 3x 


UJ + 

J 


- 2 


W ['\.E j K.e} 


(4.20) 
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After interchanging the W and h matrices equation 1+.20 becomes 


''ij = hJ [“J ^ ^ L^J H 

M{%J 


As before {(f).} is replaced by [IPHX] {<f),}. 

J J 


Similarly: 


h }= 


[EP] ((f) }, 
J 


where [EP] extrapolates (cf).} to the trailing edge. 

J 


ih.2l) 


(4.22) 


is formed from with the help of a differentiating matrix [blD^ : 




, [DHX] 


(4.23) 


h. i 

_ ^T.E._ 


is found by extrapolation similar to equation 4.22. 




Combining equations 4.21 through 4.24 the expression 


“2ik [IPHX] +2 [DHX]'^ [IPHX] 


- 2 [EH]' 


|^Wj.j [EP] {^. 


is obtained. 


(4.24) 


(4.25) 
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Letting 


[AIC] = 


-2ik 


[IPHX] +2 [DHX]' 


[s] 


[IPHX] + 


(4.26) 


[EM] 


T 


b] 


[EP] 


[VPIC] 


and following the procedure that led to equation 4.19* again an equation 
is obtained that can be written as 


[A^j] = [2]^ [H]'^ [AIC] [W] [7] 


(4.27) 


4.1.4 Conclusion., It has been shown that rather simple matrix equations 
relate the generalized aerodynamic forces to velocity potential influence 
coefficients. Wliether equation 4.7 is followed or equation 4.8, the same 
expression for the matrix of generalized aerodynamic influence coefficients is 
found when using aerodynamics based on the velocity potential as is found 
when using the integral kernel function approach (Equation 5 .12 of 
Reference l). 


4.2 Hunting Due to Piecewise Interpolation 

It is .generally accepted that when the matrix of generalized aero- 
dynamic force coefficients, [A(k)], is determined for a discrete set of 
values, kjj^ , of the reduced frequency, interpolation is adequate for 
approximating [A(k)] at arbitrary values of k. One method of interpolation, 
used successfully at the Lockheed California Company, is the cubic 
polynomial . 

For the cubic polynomial Lagrange’s interpolation formula is considered 
to be most efficient since it expresses [A(k)] directly in terms of its 

value ^A(k^) 

4 

j^A(k) [Adcj)] 

1=1 


J , at discrete values ^ 3, 4: 
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where is defined such that: 


(k-k^) (k-k^) (k-kj^) 

" Tki-k^) (k3_-k3) (k^-kj^) 


(J+.29) 


Cyclic substitution leads to and SB^. 

The interpolation formula U.28 is considered most accurate for the 
interval k^ < k k^. For the interval k^ < k < kj^ the index S, must be 
increased by one; for k^ < k < k^ the index must be lowered by one. 

Reference 1 indicates that direct differentiation of equation 4.28 leads 


to jumps in the derivative of [A(k)] at all values k = k^ if the cubic 

polynomial formula is reindexed as k moves from interval to interval. The 
j;mips in the derivative can be avoided if the differentiation precedes the 

interpolation. This is accomplished by evaluating [A*(k)] at k^ and 

determining [A'(k)] - for arbitrary k from [A'(kj^)] in the same way as 

[A(k)l is determined from [A(k^)]. These jumps in the derivatives of the 

aerodynamics matrix could lead to oscillatory non-convergence ("hunting”) of 
the structiiral resizing if the resizing procedure used requires the deriva- 
tives of the flutter speed. The problem does not occur if one interpolation 
formula is used for the entire k range or if special provisions for con- 
tinuity are made at the initial set of discrete k values: k^^. The latter 
is the case if cubic spline interpolation is used. 

A possibility of hunting exists that is not related to jumps in the value 
of the derivative. It has occurred while solving the flutter equation accord- 
ing to the p-k method (Reference 3). In the p-k method a modified one- 
dimensional Regula Falsi is used to solve for the stability roots p - (y+i) k 
at constant velocity (see also Section 2.2). In the Lockheed-California 
Company's p-k method computer program, cubic polynomial interpolation is 
used for the aerodynamic matrix. As the value of k moves from one interval 
to the next during the Iterative solution process, two different aerodynamic 
polynomials are actually used in solving the p. This can lead to hunting, 
in which the solution is sent back and foi’th between two adjacent k inter- 
vals. This type of hunting has been overcome by allowing k to take values 
that lie some distance into an adjacent interval without changing the inter- 
polation formula. Thus the interpolation formula for the i^^’tei'val 


Only when k moves out of this 


is actually used if 


e < k < + e. 


modified interval is the interpolation formula for the modified adjacent 
interval used. Thus for a k range equal to 2e, straddling each kjj^, the 
aerodynamics matrix is not uniquely defined. Actual experience with the 
modified interval is restricted to e equal to half the adjacent interval. 
This has not led to practical difficulties. 


4.3 Comparison of Options for Forming 
the Aerodynamics Matrix 

The matrix of generalized aerodynamic force coefficients was formulated 
in equation 4.19 as: 

[A(k)] = [z]^ [H]'^ [AlC(k)] [W(k)] [z] (4.31) 


where [AlC(k)], a function of the reduced frequency k = ^ and the Mach 

number, is the core of the aerodynamics. It is independent of mode shape and 
can be evaluated by any suitable aerodynamics method. Its elements are basic 

aerodynamic influence coefficients defining lumped aerodynamic forces {Z } 

a 

at an aerodynamic force grid in terms of the angles of attack at downwash 
collocation points: 

{Z } = [AlC(k)] {a} . 

The matrix. [W(k ) ] = ^[DX] + ik[DZ]J relates the angles of attack to the 

structural displacements {z}. It is independent of mode shape. 

T 

The matrix [H] is independent of k and of mode shape, and distrib- 
utes lumped aerodynamic forces and moments to the structural grid. 

The matrix [z"] is the matrix of modal columns in terms of the struc- 
tural displacements {z}. 

T 

The matrices [AIC(k)l, [H] and tw(k)] are constant during an opti- 
mization procedure. They will be used many times during the design process 
of an airplane with a given external configuration. It is therefore advan- 
tageous to form these matrices in a special aerodynamics computer program. 




4i 


(r 


Each time during an optimization procedure that a reniodalization takes 
place [A{k)l must be recomputed. Depending on the dimensions of the 
matrices in equation U.31 it may be more efficient to compute the triple 
matrix product 

[HAW(k)] = [H]^ [AlC(k)] [W(k)] (U.33) 

in the aerodynamics program, or to perform one or both of the multiplications 
[^] [H] and [W(k)][z] in the optimization program. 


In this section the possible sequences of multiplication in equa- 
tion i+.31 in the form given, and with [W(k)] replaced by [DX] + ik [DZ] , 


are examined and compared. They are identified as Options Hla, Hlb, H2a. 
H2b, H3, HUa, H^b and H5 and they are defined in Section U.3.2. In all 
options it is assumed that the derivative of the ^aerodynamics matrix with 
respect to k is needed. 

Since the aerodynamics matrix and its derivative need to be evaluated 
at arbitrary values of k, interpolation is necessary. Two interpolation 
options are considered: cubic polynomial and cubic spline. 


For cubic polynomial interpolation the value of the aerodynamics matrix 
[A(k^)] and its derivative [A’(kj^)], evaluated at discrete values k = k^^ 
are used in a Lagrangian interpolation formula: 


U 

[A(k)l = ^ (k) U(kjj)] 

£.=1 

and 

)4 

[A’(k)] = ^ [A’(kj^)] 

S,=l 


( 4 . 34 ) 


( 4 . 35 ) 


where : 


(k-kg) (k-k^) (k-k^) 

'^ 1 *'^^ “ 


( 4 . 36 ') 


and cyclic substitution leads to SE ^ and Si? ^ 
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are: 


For cutic spline interpolation the expressions for [A(k)] and [A'(k)] 


[aOc)] = [aJ + [aJ (k-k^) * [aJ jk-k^)= + [A 3 ] (k-k^)3 (4.3T) 


and 


j^A'(k)J ^^ij ^ 3 


(4.38) 


For option H3, where [AlC(k)] and [W(k)] are input as the product 
[AW(k)J = [AIC(k] [W(k)l and option H5, where [AlC(k)] and [W(k)] 

are input as the product [HAW(k)] = [H]^ [AlC(k)] [W(k)], the A quantities 
are obtained from corresponding AW and HAW quantities which in turn are 
obtained from all AW(k^) and HAW(kj^) according to cubic spline theory. 

For options Hlb, H2b and H4b, where [DX] and [DZ] are input, rather 
than [W(k)]j the formulas equivalent to equations 4.37 and 4.38 are more 
complicated as shown in Section 4.3.2. 


The matrices [A(k^)] and [A'(kj^)] in equations 4.34 and 4.35j the 
matrices po] . M and ^^ 3 ^ equations 4.37 and 4.38 and the 


A matrices used in Options Hlb, H2b and H4b are formed in the flutter opti- 
mization modilLe by appropriate pre- and postmultiplications, depending on 
the option. Always, however there is an aerodynamics input recognizable 
that is related to the set of discrete k values k^ and that is generated 
outside the flutter optimization module. 


4 . 3.1 Basis for Comparison. The count formulas in Section 4.3.2 are 
developed according to the following rules. Accuracy is not used as a basis 
for comparison. For given dimensions of the matrices used to form the gen- 
eralized aerodynamic force coefficients matrix the accuracy of the aero- 
dynamics depends on the theory used, including the method for integrating or 
lumping the pressures and interpolating and differentiating the structural 
displacements. 


1. Only the computational requirements for repetitive calculations 
inside the , flutter optimization module are determined. Preparing 
the input for the flutter optimization module is non-rep^-t'itive and 
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its impact on the total optimization is small. Note, for example, 
that interpolation hy cubic spline requires considerably more com- 
putations outside the flutter optimization module than does inter- 
polation by cubic polynomial. 

2. The number of computational operations associated with matrix-matrix 
multiplications and scalar-matrix multiplications are counted. In 
general one computational operation is one multiplication followed 
by one addition. Since an addition takes only about l/lO the time 
of a multiplication, in the few places where only multiplication 
occurs, the multiplication is counted as one operation. Additions 
are not counted when occurring alone. Generating one element of a 
real matrix product requires N operations if N is the number of 
columns in the premultiplier. Generating one complex element of the 
product of a real and a complex matrix takes 2N operations. One 
complex element of the product of two complex mati'ices takes 4 n 
operations. 

3. Matrices [H] and [W(k)] , both with . K columns, may be sparsely 
populated. This enters into the computational requirements by 
indicating that each row in these matrices has (j)N non-zero ele- 
ments. Multiplication by zero is assimied to cost no computer time. 

h. The number of computational operations (Oper) includes only the 

formation of the matrices that enter into the interpolation formulas, 
but not the operations for the actual interpolation. 

5- The number of matrix element words (one word for real, two words 
for complex element) to be stored (stor) for easy access by the 
flutter optimization module is determined. 

6, Since it may be desirable to keep the matrices involved in the 
interpolation of the aerodynamics in core during the optimization 
procedure, the required core (Cor) space in terms of the number of 
matrix element words is determined. 

7. During an optimization the value of k may drift from one interval 
to the next. If the matrices involved in the interpolation are kept 
in core, then it is of interest to know how much of the core must 

be replaced as k drifts to the next interval. This is called the 
read-in count (Read). 

h.3»2 Count Formulas. Consistent with the preceding general discussion, 
specific count formulas for several options are developed in this section. 

The following may assist in interpreting the formulas. 




The matrices [A(kj2^)] and [A*(kj^)] are the matrix of generalized 

aerodynamic force coefficients and its derj^vative to be evaluated outside the 

flutter optimization program for discrete values k = kg^. These matrices are 

to be used in a Lagrangian interpolation formula corresponding to a cubic 

polynomial. The number of k^ values required is = n^ + 3 , where 

is the number of k intervals needed. The matrices [A(k)] and [A'(k)] are 
the interpolated aerodynamics matrix and its derivative. These symbols are 
used in connection with the cubic spline options. 


The barred quantities : 


AIC,.l , AW.l , HA. 

L L Jj L -^J 



to be used 


in equations of the form given in equation 1+. 3T, are functions of all k| 
values and a set of four (j=0,l,2,3) is computed for each interval, as 
required, with the help of cubic spline formulas. 


Additional shortened notation is used: 

[HA(k)] = [H]^ [AlC(k)] , [AW(k)] = [AIC(k)][W(k)] (4.39) 

Multiplications contributing to the number of computational operations 
to be performed in the flutter optimization program are indicated by a heavy 
dot. The sequence of multiplication is indicated by Arabic numerals. Prod- 
ucts, already evaluated, but re-used in another part of the same option lack 
the multiplication dot and are assembled within matrix brackets; e.g. 

[[z]'' [H]^j . 

The dimensions of the matrices to be multiplied are suiranarized below. 

The symbol i{)ll is used when the actual dimension is N, but there is a pos- 
sibility that the matrix is sparse. 

The interpretation of the matrix dimensions used is as follows: f 

M: number of modes used in flutter analysis 

N: number of discrete degrees of freedoms: structural displacements 

K: number of lumped aerodynamic forces 

D: number of downwash collocation points 

(}): fractioh of non zero elements in rows of [W(k)] and [H]. 
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Options HI and H2 

[H]^ [AIC(k] ,[W(k)] [z] 

(M,N) («j)N,K) (K,D) (D,(t>N) (N,M) 

Optidn H3 

[zf [Hl’^ [AW(k)l ra 
■ (M,N) (P,K) (K,K) (N,M) 

Option H4 , 

[i-f [HA(k)l [W(k)] [7] 

(M,N) (H,D) (D.p) (N,M) 

Opt ion H5 

[zf [HAW(k)] [¥] 

(M,N) (N,N) (N,M) 


Option Hla CuTaic Polynomial 


[A{k^)] = [z]'' . [h]'' . [AlO(k^)] • [w(k^)] - [*] 


1 


V y. — ' 


•V 

3 


V 


|^A»(kj^)j = |^[zf[H]^j • ^AICMkj^)j- |^[W(k^)] [z]j 


"V" 

5 


"V 

6 


[ 


+ 1 [zf [nf [Aic(k^)]j 


■V 

8 


he 


ih.ko) 


ih.hl) 


(U.I+2) 


ih.h3) 


• i [DZ] [z] 



Computational operations j 


sequence 

1 

2 

3 

1* 

number 

(J)NMK 

2n^(j)NMD 

2n^KMD 


sequence 

5 

6 

7 

8 

mmiber 

2n^KMD 


([)NMD 

2n^DM^ 


Total number of operations: 


Oper Hla Poly = (jjKMK + (2n +l) ^NMD + J+n^KMD + lOn DM^ 

^ A X 

Required input quantities: [H], [AlC(k^)] , [AlC'(k^)] , [W(k^)J and [DZ] 

Stor Hla Poly = <J)NK + + (2n^+l) (f>ND 

Required for interpolation: [A(k^)] and [A'(k^)] for four k values 

Cor Hla Poly = l6 

Required read-in for new Interval: [A(k^)] and [A'(k^)J for one k value 

Read Hla Poly = 4 

Option Hla Cubic Spline 
3 

• [H]^ . [Sc ] . [W(k)l . [i-] (k-kj*^ 


[A(k)] = 





[A'{k)]=2^ I^h] [^AICjJ I^W(k)] [zj J (k-t 

t) ^ - 


[^r H" 1 • [”"] • H " 

.1=0 L J _ _ 


In the Hla sequence [W(k)] is computed for each interpolated k 
value, which is accomplished by evaluating [W(k)J = [DX] + ik [DZ] . There- 
fore. Option Hla Cubic Spline is identical to Option Hlb Cubic Spline. 

Option Hlb Cubic Polynomial 
[A(k^)] = • [H]^ • [AlC(kj^)] • [DX] • [z] + 


k^ [z]'^ [H]'^ [AlC(k^)] i [DZ] [z] j 
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A'(k^) = [z]^ [H]^ . AlC’(k^) • ^ [DX] [z] + 



Computational operations: 


sequence 

number 

1 

(>MNK 

2 

(|)MND 

3 

2n^KMD 

4 

2n,DM^ 

A 

5 

2n^KMD 

6 

2n.DM^ 

A 

sequence 

number 

7 

(J.NMD 

,8 

0 

9 

2n,DM^ 

A 

10 

2n/ 

11 

2n^DM^ 

12 


Total number of operations: 

Oper Hlb Poly = cj)NMK + 2<j)NMD + 4n.KMD + 8n,DM^ + l*n,M^ 

A A A 

Required input quantities: [H] , [AlC(k^)] , [AlC’(k^)] , [DX] and [DZ] 

Stor Hlb Poly = 4>NK + 4n^KD + PtfiND 

Required for interpolation: [A(k^)] and [A'(k^)] for four k values 

Cor Hlb Poly = l6 

Required read-in for new interval: [A(k^)] and [A’(kj^)] for one k value 

Read Hlb Poly = 
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i 
















Option HI~b Cubic Spline 


[A(k)] = ^ • [H]*^ • [AIC^] • [DX] • [z] (k-k)*^ + 

j =0 ' T " ' o ^ 


' Y ' 

3 


T 


^ |"ra'^[H]’^[AlC.]] • [DZ] • m ik (K-kj)'’ 
J=0 ■ ^ ' 


■V" 

6 


o 

[A»(k)l = ^ J^[7]'^[H]'^[nc^] [DX] [7] j j (k-k^)*^ 
j =0 ^ 


-1 


+ ^ j^[z]^[H]^[AIC^] [DZ] [z] 
J =0 


] 1 


+ '2 [ra’’[H]''[5IUj] [DZ] [?] 

3-0 


] "‘'J i^-h )■’ 


-1 


Computational operat ions ; 


sequence 

1 

2 

3 

4 

5 

6 

number 

4>NMK 

PDM 

8n^KDM 

0 

8n.DM^ 

D 

<{)MDM 

8n^DM^ 

0 


Total number of operations: 


Oper Hlb Spline = (|>NMK + 2(|)NDM + Bn.KDM + l6n.DM" 

0 0 
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Required input quantities ( for .definition of barred quantities see the 
beginning of this section and equation U.3T): 

[H] , [ATC.] , [Ainj , [Al^,] , [AT^„] , [DX] and [DZ] 
u ± i 

Stor Hlb Spline = <t>NK + 8n.KD + 2^ND' , 

o 6 

Required for :i-nterpolation : 

[F]'^[H]'^[aTC.] [DX] [z-] and [^]'^[H]'^[AIC . ] [DZ] [ z ] , 

J J 

four of each per interval (j=0,l,2,3). 

Cor Hlb Spline = l6 

Required read-in for new interval: Replace complete core 

Read Hlb Spline = l6 


Option H2a Cubic Polynomial 


[A(k^)] = [zf • [H]'^ • [AlC(k^)] • [W(k^)] • [z-] 

' y ' ^ ' 


■V 

3 


■V" 

4 


[A’(k^)] 


= [^[zf[H]'^j-j^AIC'(k^)j«|^[W(k^)][z]j + j^{z]'^[H]'^j*^AIC(kj^)j • i[DZ][z] 


T" 


■V 

6 
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Computational operations: 

12 3 

<|>NMK 2n.<}>NMD im.KDM 

A A 

5 6 7 , 

Un^KDM 2n^KM^ (j>NMD 

Total niimber of operations: 

Oper H2a Poly = (j)NMK + (2n +1) (|)NMD + lOn.KDM + 6n.KM^ 

A A A 

Required input quantities: [H] , [AlC(k^)] , tAIC'(k^)] , [W(k^)] and [DZ] 

Stor H2a Poly = (|)NK + Un KD + (2n.+l) (|)ND 

A A 

Required for interpolation: [A(k^)] and [A’(k^)] for four k values 

Cor H2a Poly = l6 

Required read-in for new interval: [A(k )] and [A'(k )] for one k value 

X X 

Read H2a Poly = 4 


2n^KM 

A 


2n^KDM 


2n^m 


sequence 

number 

sequence 

number 


Option H2a Cubic Spline 


[A(k)] = y [z]"^ • [Hl'^ • [AIC.] • [W(k)]_ • [z] (k-kj*^ 

^ t J .1 I j Jc 












[A'(k)] 


■ 2 [ 

j=0 


[z]'^[H]^[MCj] [W.(k)] [z] j 


j-1 


^ |'[i'l'’^[H]’^'] • [AlCj] • pz] •' [F]_ i (k-kj^)' 
J=0 ' > ' 

' , / 

6 


In the K2a sequence [W(k)] is computed for each interpolated k 
value which is accomplished by evaluating [W(k)] = [DX] + ik [DZ]. There- 
fore Option H2a Cubic Spline is identi'^^al to Option H2b Cubic Spline. 

Option H2b Cubic Polynomial 


[A(k^] 



T 


+ k 

V. 


.[ 


[z]^[H]'^[AIC(kj^)] i [DZ] [z]j 

Vr- 

11 
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Required input quantities: 


[H], [AlC(k^)], [AlC(k^)], [DX] and [DZ] 

Stor H2b Poly = pK + Itn.KD + 241 ND 

A 

Required for interpolation: [A(k^)] and [A'(k^)] for four k values 

Cor H2b Poly ='i6 

Required read-in for new interval: [A(k^ ) ] and [A* (k^ ) ] for one k value 

Read H2b Poly = 4 


Option H2b Cubic Spline 


[A(k)] = y [zf • [H]"^ . [AIC ] • [DX] . [z] (k-kj*^ + 

J )L 

j=0 V / 




■y" 

3 


"V 

k 


j=0 L 


[zf [H]'^ 


[AIC ] • [DZ] • [z] ik (k-kjj 
<j ■ V / I 


~V" 

6 
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[A'(k)l = 


3 r 

s 

j=o L 


[zf [H]'^ [AICj] [DX] [z] 


j + 


3 r 

0=0 l- 

3 r 

3=0 L 


[z]^ [H]"^ [AICj] [DZ] [z] 


[z]"^ [H]^ [AIC ] [DZ] [z] 
J 


i (k-k^)*^ + 


ik j (k-kjj^) 


3-1 


Computational operations: 


sequence 

1 

2 

3 

4 

5 

• 6 

7 

number 

(|)NMK 

(()NDM 

8 n^KDM 

0 

8 n^KM^ 

Q 

4 .NMD 

8 n.KDM 

0 

8 n^KM^ 

6 


Total number of operations 

2 

Oper H2b Spline = <j)IJMK + 2<j)NMD + l6ngKDM + l6ngKM 
Required input quantities: [H], [AICp^], [AIC^l, [AIC^l, [AIC^] , [DX] and [DZ] 

St or H2b Spline = (f)NK + 8n.KD + 2 ND 

0 9 

Required for interpolation: [z]'^ [H]*^ [AIC.] [DX] [z] and 

t3 

[z]"^ [H]"^ [aTc . ] [DZ] [z], four of each per interval. 

t) 

Core H2b Spline = l6 

Required read-infor new interval: Replace complete core 

Read H2b Spline = l6 











Option H3 Cubic Polynomial 



[A(k^)] 


[zf * [H]^ • [AW(k^)] • [z] 


■V 

1 




■V 

2 


■V 

3 


[A'(k^)] 


[z]^ [H]^ 


[AW'(kj^)] • [z] 

' 




Computational operations: 


sequence 

1 

2 

3 

4 

5 

number 

<j)EMK 

2n.KKM 

A 

2n.KM^ 

A 

2n-KNM 

A 

2n^KM^ 


Total number of operations: 

Oper H3 Poly = (j)NMK + 4n.KNM + Un.KM^ 

A A 


Required input quantities: [Hj^ [AW(k )] and [AW(k )] 

y, 

Stor H3 Poly = <J)NK + Un.KN 

A 

Required for interpolation: [A(k )] and [A'(k )] for four k values 

X* X# 

Cor H3 Poly = l6 

Required fead-in for new interval: [A(k^)j and [AVCk^^)] for one k value 

Read H3 Poly = 4 ■ 


5T 



Option H3 Cutic Spline 


[A(k)l 


= ^ [z]"^ • [H]'^ • [AWj • [z] 


— n ^ 


j=0 


J V. 


■V 

2 


•V 

3 


[A'(k)] = 


3 r 

3" 

jt-oL 


[zf [nf [AW ] [z] 

J 


j (k-k^)'^“^ 


Computational operations : 


sequence 

1 

2 

3 

number 

(|)KMK 

Sn^KKM 

0 

8n^KM^ 

0 


Total number of operations: 

2 

Oper H3 Spline = (}iNMK + Sn^KUM + Sn^KM 
Required input quantities: [H] , [AW^] , [AW^] , [AW^] and [AW^] 

Stor H3 Spline = (jiNM + 8n^KN 

ip p .... ^ 

Required for interpolation: four [z] [H] [AW ] [z] 

J 

Cor H3 Spline = 8 

Required read-in for new interval: Replace complete core 

Read H3 Spline = 8 
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sequence 

1 

2 

nvimber 

2n.MRD 

A 

2n^(|)NMD 

A 

sequence 

5 

6 

number 

4n^^DM^ 

(j)NDM 


Un-DM 

A 



2n^MND 


2n^^DM 


Total number of operations: 

2 

Oper H4a Poly = Un^NMD + (2n^^ + l) <J>UMD + lOn^^DM 

Required input quantities: [HA(k^)], [HA'(kj^)], [W(kj^)] and [DZ] 

Stor HUa Poly = Un^ND + (2n^+l). 4>WD 

Required for interpolation: [A(kjj^)] and [A*(k^)l for four k values 

2 

Cor H4a Poly = l6 M 

Required read-in for new interval: [A(kj^)] and [A'(kj^)] for one k value 


Read H4a Poly = 4 




















uption HUa Cubic Spline 


(A(k)] = ^ [z]"^ ■ [HAj] • [W(k>] ■ [zl (k-kp'' 

j=0 ' ' ' ' 


■V 

1 


2 


*V 

3 


[A'(k)] = 


3 

V 
z. 

3=0 L 


[z]"^ [HA ] [W(k)] [z] 
J 


j (k-k^)'^"^ + 


3 r 

j-0 L 


[?]’’ [HAj 

J 


[DZ] • [z] i (k-k^)‘ 


“Y“ 

5 


In the H4a sequence [W(k)l is computed for each interpolated k 
value, which is accomplished by evaluating [w(k)] = [DX] + ik [DZ]. There- 
fore Option Hi+a Cubic Spline is identical to Option Hi|b Cubic Spline. 


[A(k^)] 


Option H^b Cubic Polynomial 


[zf • [HA(k^)] [DX] [zj + 

V ^ ^ 


■V 

3 


kj- [ [zf 


[HA(k^)] i [DZ] [z] 


"V 

7 


] 


6o 


Ji 


[A'(k^)] = [z]^ • [HA'(k^)] • [DX] [z] + k^ • [z]'^ [HA'(k^)] 


+ [zf [HA(k^)] 

— «, 


Computational operations: 


sequence 
nvun'ber 

sequence 8 ; 

2 

number 2n^DM 2n 

A 

Total number of operations: 


1- 

2 

3 , 

4 

5 

2n-MND 

A 

<)>NMD 

2n^DM^ 

(|)NMD 

0 

8 

9 

10 

11 


2n^DM^ 

2n-M^ 

A 

2n-MND 

A 

2n^DM^ 



Oper Hl|b Poly = 2(|)NlyiD + 4n^NMD + 8n^DM"^ + 4n.M^ 

A A A 

Required input quantities: [HA(k^)], [HA'(k^)], [DX] and [DZ] 

St or H4b Poly = Un ND + 2<j)ND 

A 

Required for interpolation: [A(k^)] and [A’(k^)] for four k values 

Cor H4b Poly = l6 

Required read-in for new interval: [A(k )] and [A’(k )] for one k value 

)6 Xj 


Read H4b Poly = 4 M‘" 






Option H^b Cubic Spline 


[A('k)] = ^ [zf • [HA^l • [DX] • [z] + 

j=0 ' V- ' ' Y — 

1 2 

V — ^ 


3 r 


j=0 


[z]^ [HA ] 
J 


[DZ] • [z] ik (k-k^)* 


■V 

5 


[A'(k)] = 


3 r 

j=0 L 


[z]”^ [HA ] [DX] [z] 

cl 


j (k-k^)*^"^ + 


3 r 

0=0 L 


[z]^ [HA ] [DZ] [z] 
J 


i (k-k^ r + 


3 r 

J=0 L 


[zf [HA.l [D2] [z] 
J 


ik j 


0-1 


Computational operations 


sequence 

1 

2 

3 

h 

5 

number 

SngNMD 

4>NMD 

8n^DM^ 

0 

(j)NMD 

8n.DM^ 

O 


Total number of operations: 


Oper H4b Spline = 24>NMD + Sn^NMD + l6n^ 


DM 
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Required input quantities: [HAq], [M^], [HA^], [HA^], [DX] and [DZ] 

Stor H4b Spline = 8n,ND + 2d)ND 

0 

Required for interpolation: [zf [HA ] [DX] [z] and [z]^ [HA.] [DZ] [z] 

t) 0 

Cor H4b Spline = l6 

Required read—in for new interval: Replace complete core 

Read Hllb Spline = l6 


Option H3 Cubic Polynomial 


[A(k )] = [z]^ • [HAW(kJ] • [z] 

^ V X. 7 


[A’(k^)] = [z]^ 


[HAW'(k^)] 


Computational operations: 


sequence 


number 



2n^HM 


2n^U^M 


2n^WM 


Total number of operations: 

Oper H5 Poly = l+n.N^M + 

A A 

Required input quantities: [HAW(k )] and [HAW'(k )] 

~ Si 

Stor H5 Poly = 

A 

Required for interpolation: [A(k^)] ,and [A'(kj^)] for four k values 


Cor H5 Poly = l6 
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Required read in for new interval: [A(k^)] for one k value 

Read H5 Poly = U 

O ption CuMc Spline 

3 

[A(k)] = ^ [zf • [HAWj • [i] (k-kj^)*^ 

0=0 y ' 


3 r 1 

[A’(k)] = ^ [zj^ [&] [z] 0 

i=0 L 


Computational operations 


sequence 

1 /'■ 

. ■ . 1 2 

nxxmber 

8n-N^M i 

o 

8n^WM'^ 

0 


Total number of operad^ions: 

2 2 

Oper H5 Spline = 8n^U M + Sn^NM 
Required input quantities: [HAW^], [HAW^] , [HAWg] and [HAW^] 

Stor H5 Spline = 

Required for interpolation: [^]'^ [HAW.] [z] 

J 

Cor H5 Spline = 8 

Required read-in for new interval: Replace complete core 

Read H5 Spline = 8 



6k 








H.3.3 Comparisons . The formulas developed in the preceding section are 
tabulated in Tables U-1, i+-2 and U— 3. A summary of the comparisons is 
presented in Reference 1 and is repeated below. 

Inpu t Storage Requirements - If the number of k intervals to be used 
is three or more, cubic polynomial interpolation for arbitrary values of 
k requires less input storage than cubic spline interpolation for all 
options HI through H5. 

Core Space - For options H5 and H3 cubic polynomial interpolation for 
arbitrary k requires two times as much core space as cubic spline inter- 
polation. For all other options both methods of interpolation require the 
same core space. 

Read-In - Cubic polynomial interpolation for arbitrary k requires less, 
read-in than cubic spline interpolation as the value of k moves into an 
adjacent interval. 

Humber of Computational Operations - Cubic polynomial interpolation 
for arbitrary value of k requires fewer computational operations than cubic 
spline interpolation under the following conditions: 

options H3 and H5: if the number of k intervals is more' than three. 

options HI and HU: if the number of k intervals is more than four. 

option H2: if the number of k intervals is more than five. 

Which of the options HI through H5 is most efficient is strongly affected 
by the dimensions of the matrices. From equation U. UO it can be seen that 
if K and D are small compared with N it becomes advantageous to perform 

the multiplications [zJ^Lh]"^ and [W(k)] [z] in the flutter optimization 

module. If K and D are equal to N or larger, then it becomes 

advantageous to form the product [H]^ [AlC(k)] [W(k)l outside the flutter 
module. The relationships defining when one option is bet\ter than another 
are complicated and no simple criteria have become apparent,. 

\ 

\ 


\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 
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Number of Computational Operations 


Option 


HI a Polynomial 


Hlb Polynomial 


Hla Spline 
Hlb Spline 


ine I 
ine j 


H2a Polynomial 


H2b Polynomial 


H2a Spline 
H2b Spline 


ine I 
ine j 


H3 Polynomial 


H3 Spline 


HUa Polynomial 


HUb Polynomial 


HUa Spl 
H4b Spline 


ine 1 
ine j 


H5 Polynomial 


H5 Spline 


Number of Operations 

4)NMK + (2ng+T) (f)NMD + U(n^+3)KMD + 10(n^+3)DM^ 

PMK + 2(j)NMD + 4(n^+3)KMD + 8(ng+3)DM^ + 4(n^+3)M^ 

d)NMK + 2(t)NMD + Sn.KME) + l6n.DM^ 

O 0 

({)NMK + (2ng+T) <{>NMD + 10(n^+3)KMD + 6(ng+3)KM^ 

<J)NMK + 2(|)NMD + 8(n^+ 3)KMD + 8(n.+3)KM^ + U(n.+3)M^ 
0 0 0 

(|)NMK + 2(J)NMD + l6n^KMD + l6n^KM^ 

(|)NMK + 4(n.+3)KNM + 4(n„+3)KM^ 

□ 0 ‘ 

(j>NMK + 8n^KNM + 8n.KM^ 

0 u 

4(n^+3);^>liJ t (2n^+T) (j)NMD + 10(n^+3)DM^ 

2({)Hm) + 4(n,+3)NMD + 8(n„+3)DM^ + U(n.+3)M^ 

6 0 0 

2d)NMD + 8n„NMD + l6n^DM^ 

6 0 

4(n„+3)N^M + 4(n.+3)NM^ 

0 o 

8n„N^M + 8n.NM^ 

0 0 


M number of modes used in flutter analysis f 

N number of discrete degrees of freedom: structural displacements 
K number lumped aerodynamic forces 

D number of downwash collocation points 

(1 fraction of non-zero elements in rovs of [W(k)] and [H] 

n„ niomber of k intervals to be considered 
0 

Table 4-1: Comparison of Number of Computational Operations 

in Forming the Aerodynamics Matrices 


ji 

it 

I 
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storage Required for Input Quantities 


Option 


Hla Polynomial 


Hlb Polynomial 


Hla Spline 
Hit Spline 


ine 1 
ine j 


H2a Polynomial 



<j)NK + 4(n^+3)KD + ikn^+13) <|.ND 
4)NK + 4(n +3)KD + 2(|)ND 


(J)NK + Sn^KD + 2(|>ND 

(J>NK + 4(n.-f3)KD + (2n.+T) <I>ND 
0 0 


H2b Polynomial 


H2a Spline 
H2b Spline 


ine 1 
ine J 


H3 Polynomial 


(J>NK + U(n-+3)KD + 2(j)ND 
0 


(|)NK + 8n^KD + 2(J)ND 

o 

(|)NK + 4(n,+3)KN 


H3 Spline 


H4a Polynomial 


H4b Polynomial 


(j>NK + 8n^KN 
c 


4(n^+3)ND + (2n^+T) <|)ND 
4(ng+3)ND + 2(j)ND 


H4a Spline "I 
H4b Spline J 

H5 Polynomial 


8n^KD + 2<j)ND 

4(ng+3)N^ 


H5 Spline 8n.N 

0 

number of discrete degrees of freedom: structural displacements 

ntunber of lumped aerodynamic- forces 
number of down-'i/ash collocatjion points 

fraction of non-zero elements in rows of [W(k)] and [H] 
number of k intervals to be considered 


Table 4-2: Aerodynamics Input Storage Requirements 


6 ? 








f 


■y 
■ i 

/4 


Core Space for Interpolation and Read-in for New k Interval 


Core 

Core Space 

Read-in 

Hla Polynomial 

16M^ 

4m^ 

Hlb Polynomial 

i6m^ 

4m^ 

Hla Spline *l 
Hlb Spline J 

i6m^ 

i6m^ 

H2a Polynomial 

i6m^ 

4m^ 

H2b Polynomial 

i6m^ 

4m^ 

H2a Spline "1 
H2b Spline J 

i6m^ 

i6m^ 

H3 Polynomial 

i6m^ 

4m^ 

H3 Spline 

8m^ 

8m^ 

H4a Polynomial 

i6m^ 

4m^ 

HUb Polynomial 

i6m^ 

4m^ 

HUa Spline T 
H4b Spline J 

i6m^^ 

i6m^ 

H5 Polynomial 

16M^ 

4m^ 

H5 Spline 

8m^ 

8m^ 


M mxm'ber of modes used in flutter analysis 


Table 4 - 3 : Core Space and Read-In Required for Interpolation of 

the Aerodynamics Matrices 


5. THEORETICAL CONSIDERATIONS RELATED TO 
OPTIMIZATION WITH FLUTTER CONSTRAINTS 

5.1 Introduction 

For a better understanding and appreciation of the methods of optimi- 
zation with flutter constraints discussed in Reference 1 it is helpful to 
reduce the problem of optimization to simple terms. 
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The elementary considerations presented in this section are based on 
the assumption of satisfying a minimum flutter speed requirement for several 
flight conditions and minimizing the structural weight. 

The design variables are structural sizings, such as cross-sectional 
areas of caps, skin and web thicknesses, as well as non- structural mass 
(ballast). These design variables lead to a linear relation between the 
total mass and the design variables. It is convenient to define each design 
variable by the mass it represents. Then: 

n 

M = D (5-1) 

i=l 


where M is the total mass associated with the design variables m . 

i 

Structural elements with a stiffness defined by two design variables lead 
to a nonlinear relation between M and mj^ ; for instance, if a beam cap 
is defined by a width and a height as design variables. Another example of a 
nonlinear relation between total mass and a design variable is the stiffness 
of a constant width bending element; its mass is proportional to the cubic 
root of the design variable. Adequate finite element structural repre- 
sentation is possible without such nonlinearities and thus equation 5.1 has 
sufficient generality. 


If all structural coordinates defining a structural model are retained 
in the vibration and flutter analysis, the stiffness matrix is a linear 
function of the design variables. If for reasons of practicability a number 
of structural coordinates must be eliminated this relationship may become 
nonlinear. Any nonlinear relationship between the stiffness matrix and the 
design variables is implicit in the relationship between the flutter speed 
and the design variables, V(raj^), which is nonlinear even if the stiffness is 
a linear function of the design variables. 


If there is only one flutter speed constraint the optimization problem 
can be stated as: minimize M(mi) while the condition V(mj^) = Vp is 
satisfied, where Vp is the required flutter speed. With the help of the 
Lagrangian multiplier X this problem can be formulated as n+1 equations 
with the n+1 unknowns m^ and A: 


3 

9m. 

1 


M(m, } + A {V(m. ) - V^} 

1 1 K 


V(m.) - 



(5.2) 
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Equations 5.2 reduce to 


(5-3) 


8V 

9m. 

1 


X 


and V = 

K 


The interpretation of the equations 5-3 is that at the optimum condition 
all partial derivatives of the flutter speed with respect to the design 
variables are equal. The equation is only valid for free design variables; 
i.e. variables which, at the optimum point, have not reached a minimum or 
maximum allowable value. 


In a realistic design environment equations 5.3 cannot be solved 
directly. The usual procedure is that starting with a non-optimum design a 
sequence of design changes is generated which leads to the optimum design. 

In the following sections such sequences of design changes are discussed. 
Separate sections are devoted to increasing the flutter speed to a desired 
flutter speed with a minimum mass penalty, to minimizing the mass while 
keeping the flutter speed constant, and to considering multiple flutter speed 
constraints. 


5.2 Increasing the Flutter Speed with Minimum Mass Penalty 


Let Vy be the unsatisfactory flutter speed of an initial design and Vp 
the required higher flutter speed. The problem of increasing the flutter 
speed from Vy to can be illustrated with a system with two design 

variables mp and m 2 . Figure 5-1 is obtained from a realistic numerical 
example. The structure is a simple beam representation of the wing of a 
subsonic transport. Design variable m]_ defines the incremental torsional 
stiffness over the center one third of the span of the exposed wing and m 2 
the incremental torsional stiffness of the outer one third. Contours for 
constant flutter speed and constant total mass increment over a reference mass 
are indicated. By definition mj^ = 0 and m 2 = 0 for Vy = H 50 KEAS. The 
desired flutter speed is Vj^ = 550 KEAS. Sequences of design changes in the 
direction of the Vpj = 550 KEAS curve are indicated. Four different starting 
points are examined: 0, A, B and C, and it is assumed that the starting point 

defines minimum allowable values for m^ and m 2 . Two optimum design points 
are marked in Figure 5-1. One is a free optimum, not affected by a minimum 
allowable value of either mj^ or m 2 . It is the tangential point of the 
550 KEAS contoiir and a constant total mass line. The other optimum design 
point on the 550 KEAS contour is not an extremum in the mathematical sense; 
it is a minimum mass design consistent with the constraint imposed by the 
minimum value of ra^. 

Consider two criteria for determining the direction in the m^, m 2 

dV 

design plane: direction of steepest ascent and direction of maximum • 
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5.2.1 Direction of steepest ascent. The direction of steepest ascent is a 

dm, dm 


direction, defined by the direction cosines 
dV 


1 


ds 


and (ds^ = dm^ + dm^!) 

ds • 1 (d 


for which — — is maximum. It can be shown that this maximum occurs when 
ds 


dm^ 

ds 


9V/9m, 




dm^ 

ds 


9V/9m, 




(5.h) 


and 


/ 2 

2 


/9V 

/ ( ) 


/ \ 1/ 

\ 2/ 


( 5 . 5 ) 


Equations 5*^ define a ratio between two infinitesimal increments of 

dm^ 9V/9m. 


m^^ and m^ given by: 


dm^ 9V/9m2 

The tangent to a constant flutter speed contour is defined by. 


( 5 . 6 ) 


dm. 


1 


9V/9m, 


dirig 9V/9m^ 


( 5 . 7 ) 


Thus the direction of steepest ascent, at any point, is perpendicular 
to the constant flutter speed contours. A finite design increment step in 
that direction can be defined by: 



dm^/dsj 

^ L= As- 

1 

V”sJ 1 

dm^/dsj 


As 


|9V/9m. 


1 




2 9V/9m, 


( 5 . 8 ) 


The column 


9V/9m^ 

9V/9m, 


is called the gradient of V with respect to 


the design variables. 
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In equation 5.8 the stepsize 


is defined by the magnitude of 


As, the distance travelled in the m ,m plane. With the help of equation 
5.5 equation 5-8 can be written as: 




pV/3m^' 


which relates the step size to an increment in flutter speed. 

Design paths following the direction of steepest ascent are indicated 
in Figure 5~1* It is noted that they do not lead to the optimum design 
point . 


5 . 2.2 Direction of maximum 


The direction of maximum “ is the 

dm 


direction in which the increment of the flutter speed per unit increment of 
total mass is maximum. 

At an arbitrary point in the m^jm^ design plane there is no direction 
'^l dV 

for which — is an extremum. Consider, for instance, point A as an 

initial design point and move along an M=constant contour towards the 

dV 

m^ axis: V increases, M remains the same, thus ^ ^ design change 

along an M=constant contour, however, is not feasible in view of the 
original assumption that at the initial design point mj_ and m 2 are at 
their minimum values. From simple geometric considerations it follows that 


the maximum feasible value of is in the direction of the coordinate for 

dM 

3V 3V 9V 

which IS largest. At point A and a design change parallel 

to the m^ axis results in the largest gain in flutter speed per unit mass 


added until the locus for which 


3V ^ 

am^ 3m^ 


is reached. Since this locus 


connects points representing minimum total mass designs as a function of 

flutter speed, following this locus leads to the optimum design at 

= 550 KEAS. 
ti 


dm 

Along this locus the direction is determined from the requirement 

. . 32v _ 92v _ dm 

2 2 




, from which ^ dm^ + . 


3m., 


1 2 


2 1 


3m^ 


and 


dm^ 3^Ar/3m2 - 3^V/3m^3m2 

^2 3^V/3m^ - 3^V/3m23m^ 


(5.10) 


Paths following maximimi feasible from the initial design points 


0 and B are also indicated in Figure 5-1. 

Bote that an initial design defined by point C, leads to a minimum weight 
design at which V = 550 that is not an extremum condition, due to the 
minimum value imposed on m^. 

5.2.3 Discussion . Figure 5-1 shows that following a direction defined by the 
gradient of the flutter speed, in the case shown, leads to a design configur- 
ation close to optimijm. This result cannot be generalized. In a multi-design 
variable system it is not unlikely that several design variables remain at 
their initial value for the otpimum design. If that is the case the result of 
following the steepest ascent direction generally will be farther from the 
optimTom design than is indicated in figure 5-1- 

Generalizing equation 5-9 to n design variables leads to the following 
expression for the column of design variable increments: 



Generalization of the path of maximum approach to many design 

variables takes the following, albeit impracticable form. 


Assume, as before, that at the initial design point all design variables 
are at their minimum allowable value. Obtain the largest gain in flutter 
speed per unit mass by increasing the design variable corresponding to the 


largest 


. As this design variable is increased it most likely becomes 


less effective in raising the flutter speed and after a certain amount of 
change there may be two design variables with the largest value of 


7 — . Then these two variables can be changed according to the relationship 
dm 


given by equation 5-10. When three 


's are equal the three associated 


variables can be changed in a ratio that follows from a reasoning similar to 

3V 

the reasoning that led to equation 5.10. Successively more -; 7 — 's 


become equal. When the desired flutter speed is reached the 


's associated 


with the design variables that have been increased above their minimum value 
are all equal, as was found before as a condition of optimality 
(equations 5-3). 

It should be remarked that the procedure is only presented to provide 
some insight in the process that leads to an optimum design. 


Although it is noted that a path of steepest ascent does not lead to 
the optimum design, it does lead to a favorable initial configuration from 
which to start minimizing the total mass while keeping the flutter speed 
at Vj^ . It was found (Reference 1, Appendix A) that an efficient approach 
to defining such an initial configuration is to use a resizing column 


{Am^} = C 




(5.12) 


where C is determined by requiring that the flutter speed is raised to 
in one step. It leads to an initial weight only little higher than that 
which results from a multi-step approach, but at a considerable reduction in 
computational effort. 


The presentation in the preceding sections is considerably simplified by 
defining the design variables in terms of mass units. In Reference 1 


it is shown tha,t the distribution in a resizing coliomn defined by 


3P. ’ 

X 


where P. is a generalized design variable: P. ~ » depends on the 

^ ^ i 

individual scaling of the design variables. E.g. for skin thickness 


as a design variable is different from a for cap area as a design 

variable. Since obviously the best distribution in a resizing column should 
be independent of the scaling, the effect of scaling must be eliminated by 
a normalization. The simplest normalization is to let = 1. 


5.3 Minimizing the Total Mass at Constant Flutter Speed 


5 . 3.1 Elementary Steps . Minimizing the structural mass vhile keeping the 
flutter speed constant has considerable practical significance. In 
Section 5*2 it was shown that the rather straightforward method of steepest 
ascent does not lead to the optimum design. However, the method of following 

dV 

a path of maximum is rather intricate. 

^ dM 

There are practical methods of minimizing the structural mass which may 
be used once the flutter speed has the desired value. Several are discussed 
in detail in Reference 1. It is instructive to consider an elementary 
approach, based on simple reasoning, since it provides insight into practical 
methods . 


The elementary approach is based on the following basic design change: 
take away material where it causes the smallest decrease of the flutter speed 
per unit mass removed and add material where it causes the largest increase 
of the flutter speed per unit mass added. 

The mathematical formulation is as follows: 

Let the flutter speed, V, be given by: 

V = V(m. ,m.,,..m ) = V„ (5.13) 

X ti: n i\ 


Then 


T — dm^ + - — dm, 
9mj^ 1 Sm^ £ 


= 0 


(5.14) 


Let 9V/9m^ > 0 and BV/Sm^ > 0, 3V/3m^ the largest and 3V/9m2 

the smallest of all 3V/3m's. Reducing m^ to m^-Am^ causes the smallest 
decrease in flutter speed per unit mass removed. Prom equation 5.l4 it 


follows that the amount to be added to to maintain approximately 

constant flutter gpeed is: 


av/smg 

^™1 "" 9 V / 3 m ^ ^"^2 


(5.15) 


The total mass change is 



which is a reduction. 


(5.16) 


Small steps Am^s Am 2 can be taken in succession and all derivatiyes 
recalculated after each step. As long as 9V/3m3_ and 9V/9m2 remain 
the largest and the smallest of the 9V/9m's the same design variables are 
Qhanged. As soon as a third derivative 9V/9m3 equals either 9V/9m]_ or 
9V/9mi2 it must be considered in determining the next step in a manner 
similar to what is discussed in Section 5.2.2. If m 2 reaches its minimum 
allowable value the design variable with the next to the smallest va3.ue of 
9V/9m is used instead. It is not difficult to see that a stage is reached 
in which the 9V/9m’s associated with design variables that are not at their 
minimum allowable value are all equal. Nothing can be gained by increasing 
one and decreasing another variable by a small amount: the optimum design is 

attained. 

9V , 

Of interest is the case in which some -;r — 's are negative. With 

dm 

9V/9ra2 being the smallest 9V/9m this implies 9V/9m2 < 0 and decreasing 
m 2 increases the flutter speed. The initial resizing steps, therefore, 
should be to reduce m 2 and let the flutter speed increase. Successive steps 
of reducing mp can be taken until occurrence of either one of the following: 
m2 reaches its minimum allowable value, 9V/9ni2 equals another negative 
9V/9m, or 9V/9m2 has become positive and as a result the gain in flutter 
speed is cancelled. 

5.3.2 Composite Resizing Columns. It is recognized that the elementary steps 
described in the preceding section as well as those described in Section 5.2.2 
are an inefficient approach to optimization. It is also recognized that a 
resizing column proportional to the gradient of the flutter speed (equa- 
tion 5 . 12 ) to a certain extent satisfies the basic notion of adding most 
material where it causes the largest increase of the flutter speed per unit 
mass. However, material is also added where, according to the reasoning in 
the preceding section it should be subtracted. This is indeed the very 
reason why following the direction of steepest ascent does not lead to an 
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optimuia. What is needed is a resizing column that removes the most mass where 
it causes the smallest decrease of the flutter speed per unit mass. A logical 
candidate is a resizing column that is the sum of two linearly independent 
distributions and leads to the desired effect. Equation 5 *17 shows a simple 
approach to following this concept: 


{Am^} 


_ , f 8V 1 „ r/3V\ 3V 1 

1 J ' max X J 


This can be written as: 


{Am. } 

1 


= (A + B) 


_ B 


(5. IT) 


(5.18) 




(5.19) 


which is a form identical to the resizing column for the gradient projection 
search at constant flutter speed of Rudisill and Bhatia (Reference 6 and 
Reference 1, Section 6.2). 

Kdcs removal where it causes the smallest decrease of the flutter speed 
per unit mass is emphasized if the second column of equation 5.17 is 


replaced by 


1 _ 

' 9V/9m. ■ 


“i» = A ji^j- b|,^J 


( 5 . 20 ) 


The constant A is related to the total mass addition, W^^ , associated 
with the first column: 


A = M 


|_lj {9V/9m^} 


( 5 . 21 ) 


The constant B can be determined from the condition 




{Am^} = 0 


( 5 . 22 ) 


which would result in zero change in flutter speed if V were a linear function 
of the design variables mj^. 
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Applying equation 5-22 to equation 

n 

i=l 

Due to the presence of nonlinearities, however, equation 5.23 does, in 
general, lead to a finite change in flutter speed. 

A condition somewhat more complicated than equation 5.22 follows from 
the requirement that the flutter speed is to be kept exactly constant. It is 
determined by solving equation 5.2U for k and B according to the method 
of Incremented Flutter Analysis 


(5.2li) 


or 


The value can be chosen arbitrarily or it can be based on a step 

size criterion, such as discussed in Reference 1, Section 6.6. 


( W {3V/9m. } 
g. Y. k, V, fall. [aY/3n,^> » B fVM) 

where {VM} corresponds to the second column of either equation 5. IT 
equation 5-20. 



5.20, for example, leads to: 




(5.23) 


The resizing coliomn defined by equation 5.20 has been used successfully 

8V 

for a case in which all s had positive values. Results are presented 
in Reference 1. 


When there are negative 8V/9m's the 9V/9mi contribution in 
equations 5.17 and 5-20 tends to remove most mass where it increases the 


flutter speed most. On the other hand the column in equation 5 <20 

9V ^ 9V 

adds to design variables with negative — — and tends to outweigh the — — 

dm dm. 

X 

contribution. 


Figure 5-2 shows, graphically, possibilities of modifying the columns 
in equation 5*20 to accommodate negative 9V/9m's. The elements are 
arranged in the order of increasing value of 9V/9m on the vertical axis. 
The value of each element is measured along the horizontal axis. 


Figure 5-2 (a) shows, for reference, the components of the resizing 
column of equation 5-20, but for positive 9V/9rn’s only. Figure 5-2(b) 


extends into the range of negative 


9V 

9m 


s. For negative values of 9V/9 


m 
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/ 



value of element value of element 

= av/9m 


dV/Bm 

(a) Positive 9V/ 9m' s only- 



value of element value of element 

= 9V/9m 

(b) Replace Negative 3 ^ by Zero 



(c) Two Non-Negative Columns 



value of element value of element 

= ( - 8v/a 

(d) Avoid near Infinite Elements 

Figure 5-2: Composite Resizing Col-umns, 


the elements in the second column are set to zero to avoid the problem 

associated with large negative values of — i — . 

9V/3m 

Figure 5~2(c) is a variation of Figure 5~2(b) in that the negative 

— s have been moved, with opposite sign, from the first to the second 

column. Thus neither column has negative elements. This has advantages 
for including minimum size constraints in resizing algorithms. 

Equation 5-lT is represented in Figure 5-2(d), also with the modifica- 
tion that avoids negative elements in the component columns. 

Since in an actual airplane design several or many 3V/9m's may be 

very close to zeroj the column { 9V/9m. ;} may contain extremely large element 

If this would constitute a problem a limit could be set on the value of 

> or the problem could be avoided by using the resizing colxomn of 

Figure 5-2 (d). 


s . 


5.4 Multiple Flutter Speed Constraints 

5-4.1 Introduction. When the flutter analysis of an airplane design indi- 
cates significantly unsatisfactory flutter characteristics, it often happens 
that for different loading conditions (full fuel, no fuel, full pay load, 
etc.) and/or Mach number several instances of flutter speeds less than the 
required Vp occur. Even for one loading condition and one Mach ntimber 
several zero modal dampings at speeds below Vp may occur. When resizing the 
structure to satisfy the flutter speed constraints one may think of trying to 
move all zero damping intersections up to Vp. Whether or not this is 
possible is irrelevant from a practical point of view, since it does not 
necessarily lead to the lightest structure, as is shown in the next paragraph, 

Assiune a structure with tuo flutter speeds, 'and Yq , both below 
Vp. In some fashion (uniform increase of stiffness level, steepest ascent, 
or other) the structure is modified such that Vp = Vp. Assume that now 
Vp > Vp. A one-flutter-speed-constraint mass minimisation, keeping Vp = Vp, 
can now be performed. Consider the case for which V 2 remains above Vp 
until uhe minimum mass for Vp = Vp is reached. Then tha,t is the minimimi 
mass for which the flutter constraints Vp > Vp and Vp ^ Vp are satisfied. 
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To lower Vg to while keeping V]_ = Vr, if it is at all possihle, 

requires an increase in mass since the minimum weight for which = Vp is 
already attained. 

If in minimizing the mass for Vq = Vp, Vp decreases to Vr, then hoth 
flutter speeds should he maintained at Vp until a minimum mass condition is 
reached. 

When there are more than two flutter speeds to contend with, minimizing 
the mass is usually accomplished by the following procedure. 

Considering the complete set of flutter characteristics (complete f-g-V 
diagrams for various . Mach numbers and loading conditions) of the initial 
structure, it is estimated which of the zero damping intersections is most 
critical. The structure is modified to move this intersection to Vj^. If 
there still remains an intersection at a speed below Y-^ the structure is 
further modified. The result is an initial non-optimum structure with 
Vi = Vr and Vx>Vr (i=2,3,..v). A one-flutter-speed-constraint mass mini- 
mization, keeping Vp = VR, is performed until a minimum mass design is 
reached, or until another flutter speed, say Vp , decreases to Vp = VR. In 
the latter case the weight minimization is continued with two flutter speed 
constraints: Vi = Vp = VR. Eventually a third and more speeds may be drawn 

into the minimization procedure. 


The following section provides mathematical background for the above 
discussion. 


5.h.2 Mathematical background . If two flutter speed constraints are active 
the optimization problem can be stated as; minimize M(mi) while the 
conditions V2 (mx) = Vri and Vp(mj^) = Vr 2 are satisfied. (Note that dif- 
ferent flight conditions may have different -Vr). With the help of the 
Lagrangian multipliers and Ap this problem can be formulated as n+2 

equations with the n+2 unknowns m^, A]_ and Ap: 


8 

3m. 




= 0 


V2(m.) - = 0 


(5.25) 


These equations reduce to 


3V 3V 

= ° i=l,2,...n (5.26) 

1 X 
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and 


V = V 
1 R1 


V = V 
2 R2 


(5.27) 


Consider three of the equations 5.26, say for i = 1, 2 and 3; 


9V 8V 

-r — + T — +1=0 

1 2 9in^ 


3V 9V 

A, ^ — + A- +1 = 0 
1 2 


3V 9V 


( 5 . 28 ) 


Only non-zero values of A and A can satisfy equation 5.26 for 
9V^ 9V2^ 

finite values of and . The condition for non-zero solutions for 

i i 

A^ and A^ follows from equations 5 . 28 : 


1 

1 

1 




9m^ 

9V 

9V 




1 

1 

9m^ 

= 0 

(5.29) 

‘ti 

ii 

■av2 

av2 





9m^ 

9-2 

9m2 





Since the subscripts 1, 2 and 3 are arbitrarily assigned, equation 5.29 
must be satisfied for all possible combinations of three design variables. 


[n the case of one flutter speed constraint equation 5.29 reduces to: 


(5.30) 


from which: 


3V, 


3V. 


3m-, 


3m^ 


1 “"‘2 

which corresponds to the first equation of equations 5.3. 

Extension of equation 5-29 to v flutter speed constraints leads to: 


1 


3ra, 


C 

3m, 


1 

3m^ 

3V. 

C 

3m^ 


1 

3V„ 


3m 


v+1 

3V^ 


3m 


v+1 


3V 



3m, 


3V 

I 

3m^ 


3V 


3m 


v+1 


= 0 


(5.31) 


which must he satisfied for any set m^, mg,..,m^_^^ within m^,m 2 ,..m^ 
(n > v+1). 

Equation 5.31 is satisfied if any row of the determinant has equal 
elements . 


3V^ 3V^ 3V^ 

3m, ~ 3m^ 3m 

12 n 


(5.32) 


This is the minimum mass condition for one flutter speed constraint. _ It 
means' that if the mass is minimized while satisfying one flutter speed con- 
straint Vp = Vpp and all other flutter speeds exceed their associated Vr, the 
minimum mass for 2 Vfji (i~l»2>..v) is obtained regardless of the values 
of the derivatives in the other rows. This conclusion was reached by 
informal reasoning in Section 5.A.I. 


8it 


Another special condition for which equation 5-31 is satisfied is: 



9V„ 



9V 

9V, 


2 _ 

k ^ + k 

2 . 

- k , + k^ 


1 9m^ 2 

9m^ 

1 9m^ 2 

9m^ 

1 9m 2 

9m^ 

1 

1 


2 

n 



(5.33) 


where and k^ are arbitrary constants. 

This condition occurs when the mass is minimized while keeping two 
flutter speeds at their corresponding Vj^. If all other flutter speeds equal 

or exceed their associated , then the minimum mass for V. ^ V . 

K 1 

(i=l,2,..v) is obtained regardless of the values of the derivatives in the 
rows not represented in equation 5-33. The reasoning can be expanded to cover 
any number of antive flutter speed constraints. 

The same restriction is imposed on equation 5-31 that is imposed on 
equations 5.3: equation 5*31 is only valid for free design variables; i.e. 
design variables which, at the optimum point have not reached a minimum or 
maximum allowable value. 

5.4.3 Discussion, Considerable background on optimization with one flutter 
speed constraint has been published and more is provided in Reference 1. The 
elementary considerations related to the one-flutter-speed-constraint opti- 
mization in preceding sections, therefore, relate directly to published 
material. 

Although several methods of optimization appear sufficiently generally 
formulated to be applicable to the case of multiple flutter speed constraints 
the literature does not provide any examples. The following discussion, 
therefore, is rather speculative. 

5. 4. 3.1 Increasing the Flutter Speeds 

In the case of one flutter speed constraint a simple resizing column is 
defined, proportional to the initial gradient of the flutter speed (equa- 
tion 5.12), to increase the flutter speed to in one step. Numerical 

examples have shown that this provides a favorable initial design from which 
to start a mass minimization process at constant flutter speed. 


The resizing distribution defined by equation 5.12 follows the notion of 
adding to each variable an amount proportional to its effectiveness in 
increasing the flutter speed. Extending this notion to the case of several 
unsatisfied flutter speed constraints a resizing column could be conceived in 
which each element is proportional to the effectiveness of the assocxated 
design variable in increasing all submarginal flutter speeds. This could 
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lead to a resizing column that is proportional to the sura of the gradients of 
the flutter speeds: 


{ Am . } = C •< 
1 


u 




3V 

a. 

3m. 

1 


(V. = 


flutter speeds helow V ; j = l,2..u) 

R 


(5.34) 


However, one flutter speed may he close to its Vp, another far below. 
Obviously, design variables favoring the latter should be weighted more 
heavily than variables favoring the former. This leads to a resizing column 
based on a weighted sim of the gradients of the flutter speed: 


{Am. } = C 
1 


A / \ 9V. 

z K. -3“ 

j=l ^ J / - 


(5.35) 


The scalar C should be determined such that one flutter speed V =V 

1 R1 

and all other flutter speeds V. > V . (j=2,3..u), 

j Hj 

It is emphasized that the authors are not aware of any numerical evalu- 
ation of the suggested procedure. 

5. 4. 3. 2 Minimizing the Mass 

Consider fii’st the case of one active flutter speed constraint: 

h = V ’ ''j ^ ''rj 


Use the resizing column of eq^uation 5-20, with A replaced by the value 
given by equation 5.21, and compute the resulting increment of the flutter 


speeds V. (j^l): 


A 

3V^ 3\ 


1 1 


UJ {SA/S™!} 

3m. 3n 

L iJ 


l^”iJ 


( 5 . 36 ) 


If all AV^ > 0 the optimization can be treated as a one-flutter-speed- 
Gonsti alnt mass minimization. If there is one flutter speed, say ^2* 
which AV < 0 there is a possibility that the resizing column of 
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equation 5*20 reduces below V , and a condition for can be 

^ 1 

established by equating AV^ to its largest allowable negative value: 


V - V 
R2 2 


= V - V = ■■■ , 

2 "R2 ''2 [IJ {9V 


1^1 feJfr}- • [^Jfe} 


(5.37) 


Since B is an implicit, non-linear function of (equation 5. 2it), 
an iterative process, such as a one-dimensional Regula Falsi approach, should 
be Used to determine from equation 5.37. If satisfying equa- 

tion 5-37 is smaller than the Wq that follows from the step-size criterion 
associated with equations 5*20, 5.21 and 5-2U the former defines the value 
of Wq to be used for the first resizing. .Otherwise the resizing is based 
on the criterion used in the case of one active flutter speed constraint. 

If us follows from equation 5-37 defines the step size, by defini- 


tion V = V 

1 R1 


and . For further weight reduction, while keeping 


and constant, the following may be considered. 

In accordance with the basic resizing concept expressed in Section 5.3.2: 
add where it does the most good and subtract where it does the least harm, and 
taking equation 5-20 as a model, the following resizing column can be formed: 


{Am. } = A 
1 


3m. 

L 1 


(5.38) 


The requirements — 0 arid AV^ = 0 lead to the following two com- 

plex determinantal equations which comprise four equations with the four 
unknowns: two reduced frequencies k^ and k^ ,: and A and B. 


^"^ 2 ! f 1 l\ 

oY, ^ I ° 

3m.. ^ 3m. / 

I 1 iJ / 


( 5 . 39 ; 




g.Y,k 2 ’^R 2 ’'^ 

\ 


_ f 1 1 \ 

V”i *“ij ^ ^''2 ■ 1 

Bm . 9m . / 

I 1 1 J / 


= 0 


:5.uo) 


Although;.:equations 5-39 and 5.^0 suggest unique solutions for A and B, 
and the k values associated with and Vj ^2 > these solutions may be 

complex values and thus without physical meaning. 

It is noted that the two complex equations 5.39 and 5.40 with the four 
unknowns A, B, k^ and k 2 can only be formulated because of the non- 
linear character of the relationships between the flutter speeds and the 

design variables. Keeping AV^ = 0 and AV^ = 0 on a linear basis leads to: 


9V^ 

/ 

r 9 V av 1 

1 

/a. 

' 

CM 

H 

r— 1 j 

9m. 

9m. 9m. I 

1 


L 1 i.J 



/ 

r 9V, av,,'i 

2 

/a. 

1 + 2 , 

9m. 

9m. 9ra. 
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L 1 ij 


V- B 


9V, 



1 


2 

9m. 


9m. 

L 1 


1 

r. 




1 


9V, 



1 


2 

9 m. 


9 m. 

L 1 


1 


= 0 


( 5 . 41) 


= 0 


(5.42) 


which are two homogeneous equations with essentially one unknown: A/B, and 

in general, no solution for A and B separately. 

It is also noted that, if real values of A and B are found, the 
distribution of {Am^} as well as its magnitude is defined. This may well 
lead to convergence difficulties. 

A more promising approach is to select three linearly Independent col- 
umns, such as shown in equation 5-43 or equation 5 . 44 . 


{Am^} 


= A 


1 9m^ 9m^ j 


B 


{Am. } 
1 


= A 


9V, 


^ + 

2 

3m. 

9m. 

1 

r 1 



+ C {1} 


( 5 . 43 ) 


'dY 9V " 

+ B r 1 1 ^ r i i 

^9mi 9m J jaV^/Bm.j: 


(5.44) 
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5.4.3. 3 The Distribution 


'am.’’ 8m. ( 

L 1 ij 


Support for the use of the distribution 


as a guide for 


A valiie of A, defining stepsize, can t>e selected, and equations equiva- 
xeht"'to equations 5.4l and 5-42 make it pos| 5 ible to relate B and C to A. 
Since equations 5.4l and 5*42 are linear, and V 2 are not maintained 

exactly at 


generating a resizing column is found in the feasible directions procedure of 
Reference I 6 . It has been shown (see discussion in Reference 1, Section 6 . 5 ) 
that with the procedure suggested in Reference I 6 (push-off factor 6=1) 
the distribution of the resizing column can be determined from the following 
equations : 


smallest 


— ^ {S} - 


(5.45) 


largest 


{ft}"”} 

■ft}-4‘" 

J 


(5.46) 


(5.4t) 


The values of the elements S^^ and follow from the two equatipns 

with two unknowns: 5*46 and 5 •4?. 
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9V^ 

Thus design variables for which t — + -r — is large are increased; those' 

dm dm 

3Vi av^ 

for which -r — + — — is small are decreased, 

am 3m 

Numerical results thus obtained were in exact sgreement with results 
obtained using the Simplex method as suggested in Reference l6 


5.5 Optimality Condition, Ballasting and Stiffening 


The optimality condition for one flutter speed: 


3m^ am^ 


av ^ 

am^ am. 


:5.u8) 


for free design variables is discussed in connection with practical considera- 
tions regarding ballast and stiffness design variables. 

The literature pays little or no explicit attention to ballast (dead 
weight) as a design variable. The reason for this omission is understandable: 
any method of optimization that can handle design variables representing 
related stiffness and mass changes can handle a design variable representing 
a mass change only. 


The need to consider mass ballast as a design variable is an additional 
argument in favor of using mass as the dimension of design variables. 

It should be noted that a stiffness change with zero associated mass 
change cannot be represented by mass units. Such a stiffness change, however, 
must have some other penalty such as drag or production cost, and the objec- 
tive function to be minimized cannot be the sum of structural and ballast 


Inclusion of design variables representing ballast obviously increases 
the number of design variables. Although ballasting design variables are 
less costly in term of computer time and problem formulation than stiffness 
design variables, an indiscriminate increase is an unwelcome burden. In 
Reference 1 it is suggested that it may be convenient to first optimize the 
structdle using stiffness design variables and their associated mass only and 
then to determine whether any mass change by itself is more efficient than the 
stiffness changes in raising the flutter speed. Alternately, the values of 
flutter speed derivatives with respect to ballast mass, determined at a stiff- 
ness^ optimum, may be used in an attempt to further lower the total mass. The 
following procedure could be followed. ' 
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- l,2,...,n^) represent stiffness design variables and 
nig(3 = n^+2,...n) ballast design variables. Then, after stiffness 

optimization: 


9V _ 8V 


3V _ . . . _ 3V 

Sm 3m 

a n 

a 


(5.i^9) 


for those m^ not constrained by a sizing limit. 


Any m^ 


for which 


3V ^ 3V 

3mo ^ 3m 
p a 


(5.50) 


can not be used to reduce the total mass if it is assximed that initially there 
is no ballast present . This follows from the fact that equation 5.50 implies 
that the total mass can be reduced if mg can be reduced and m increased. 
If, however, there is a ballast configuration for which ^ 


3V ^ 3V 

3m„ 3m 

3 a 


(5.51) 


then a mass reduction by ballasting is possible. 

Thus, after an optimization in which only stiffness design variables are 
used, a new optimization may be initiated that includes all stiffness vari- 
ables and those ballasting design variables that satisfy equation 5.51. 

Although sizing limits are usually thought of as lower limits, upper 
limits can also occur. Ballasting, for instance, has an upper limit defined 
by space available and specific weight of acceptable ballasting material. 

Even structural members may have upper limit constraints; for instance in 
areas where control surface actuators must be housed in a wing box. Thus the 
optimality condition of equation 5.^8? when taking note of the thought 
expressed by equation 5.51 can be generalized as follows: 

Let m^j^ be a design variable at its lower sizing limit, m^ a design 
variable at its upper limit, and m^^ a design variable satisfying the condi- 
tion of equation 5.^8 then the generalized optimality condition is: 


3V 

3m. 

X 


= C 


3V 

3m„ 


< C 


3V 

3m 


> C 


(5.52) 


u 


The conditions 5.52 can be used to check the result of an actual opti- 
mization process. 
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6. PBACTICAL CONSIDERATIONS RELATED TO OPTIMIZATION 
WITH FLUTTER CONSTRAINTS 

Concurrent with the present study the Lockheed-California Company con- 
ducted work on contract NASl-12288 with NASA, Langley Research Center (Refer- 
ence 17 ) • Contract NASl-12288 is a design study for a supersonic transport 
with an arrow wing planform (Figure 6-1). During Task I of that contract 
three design concepts were evaluated, which included the determination of the 
weight penalty for satisfying the flutter constraints. On 'the basis of the 
results of Task I a design concept for Task II was chosen and further refined. 
One of the investigators in the present study participated in the weight mini- 
mization part of the arrow wing contract, with benefits for both studies. 

Of importance to the present study is that the arrow wing contract pro- 
vided a realistic design environment for flutter optimization. 

The approach to flutter optimization consisted of a judicious choice of 
the most critical flutter condition followed by the actual optimization by 
means of an interactive Computer Graphics program, using the method of 
Incremented Flutter Analysis to determine the sensitivity of the Flutter 
speed to design variables . 

The practical experience gained can be demonstrated by examining the 
significant findings during the Task I design of the arrow wing based on the 
ehordwise stiffened design concept. 

In the following the structural model and the flutter analyses of the 
original strength design are discussed. The structural regions delineating 
design varia,bles are defined and the non-linear relations between elements 
of the stiffness matrix and the design variables are demonstrated. Finally 
the design experience as related to the number of modes used in the flutter 
optimization and the modal updating is presented. 


6.1 Structural Model and Vibration Analysis 


The Task I structural model of the wing is symmetric about a midplane 
parallel to the x-y plane. Vertical deflections of the nodes on top and 
bottom surface of the wing are asstimed equal and in the same direction, defin- 
ing the overall deflection of the wing. Deflections in the horizontal plane 
on top and bottom surface are assumed equal and opposite and thus there are 
no elastic deflections parallel to the x-y plane in the mid-plane of the « 
wing. At each node the structural degrees of freedom are translations in the 
directions of the three axes. The fuselage is represented by a simple beam. 

The total number of structural degrees of freedom per side is 10if2. For 
symmetric modes 250 of these are translations in the z direction and, on 
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the fin, are translations in the y direction, describing the overall 
motion of the lifting surfaces. For antisymmetric ’modes 25 z translations 
on the fuselage are replaced by 25 y translations. The vibration and 
flutter equations are formed on the basis of 188 degrees of freedom for sym- 
metric analyses (Figure 6-2) and I 78 degrees of freedom for antisymmetric- 
analyses. Consistent inertia and aerodynamics matrices are based on the mops 
than 250 out-of-plare, translations. A static reduction is used to reduce the 
structural degrees of freedom from 1042 to I 88 (1T8 antisymmetric). A corre- 
sponding reduction is used for the inertia and aerodynamic matrices (Guyan 
reduction. Reference I 8 ). 

Vibration analyses are conducted to determine the lowest 50 natural fre- 
quencies and mode shapes. Table 6-1 summarizes the vibration characteristics 
of the initial, strength designed, configuration. 


6.2 Aerodynamics and Flutter Analysis 

Doublet lattice theory (Reference 9) is used to form 264th order aero- 
dynamics matrices for M = 0.6 and M = 0.9j using 233 aerodynamic boxes 
per side on the wing and 15 on the fin. They are subsequently reduced to 
188 th order by the Guyan reduction method (Reference 18 ). 

The degrees of freedom in the flutter analysis are reduced by post and 
pre multiplying the flutter matrix by a modal matrix and its transpose. 

Twenty vibration modes are used in the flutter survey of the initial 
configuration. 

Sample results of the flutter analysis are shown in Figures 6-3- through 

6 - 6 . 

The use of twenty modes was based on engineering judgment in reaching a 
compromise between accuracy and cost. Cost considerations included not only 
the direct cost per flutter analysis, but also the cost of assuring that the 
number of modes chosen, if low, is adequate. 

An example of the sensitivity of the flutter characteristics to number 
of modes used is presented in Figure 6-7- For this condition it appears 15 
modes would have given acceptable results for the flutter analysis. Regard- 
less of this answer, however , twenty modes was considered the absolute mini- 
mum nxunber of modes required /for a meaningful optimization. The effect of 
increasing the number of modes on one flutter speed is shown in Figure 6 - 8 . 

To increase the understanding of the nature of the flutter modes, modal 
participation coefficients were plotted to make visible the amount of partici- 
pation of each of the twenty natural vibration modes in in-flight modes of 
interest as a function of speed. Examination of these plots revealed that the 
flutter mode labeled bending and torsipn-v^nd the one labeled hump mode 



Figure 6 - 2 : Symmetric Degrees of Freedom for Vibration Analysis. 



MODAL FREQUENCY - HERTZ | 

MODE DESCRIPTION 

SYMMETRIC 

ANTISYMMETRIC 


OWE 

FFFP 

OWE 

FFFP 

WING 1 ST BENDING 

1.009 

0.933 

1.228 

0.908 

FUSELAGE 1 ST BENDING 

1.381 

1.206 

1.998 

1.949 

ENGINE PITCH IN PHASE 

l.6tl 

1.627 

1.514 

1.457 

ENGINE PITCH OUT OF PHASE 

1.817 

1.815 

1.821 

1.805 

FUSELA® 2 ND BENDING 

2.784 

2.261 

3.370 

3.056 

WBIG IBT TORSION 

3.288 

3.104 

3.034 

3.056 


Operating Weight Empty (OWE) - 321,000 LBS 
Full Fuel Full Payload (fFFP) - 750,000 LBS 


Table 6 - 1 : Summary of Vibration Characteristics; Initial Design. 
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FREQUENCY-, HERTZ 


SYMMETRIC FLUTTER ANALYSIS - CHORDWISE STIFFENED ARRANGEMENT 
MACH NO. = 0.6 
WEIGHT = 321 .000 LBS 




Figure 6 - 3 : Symmetric Flutter Analysis - Mach 0.6 - Operating Weight Empty. 


SYMMETRIC FLUTTER ANALYSIS - CHORDWISE STIFFENED ARRANGEMENT 
MACH NO. = 0.9 
WEIGHT ^ 750,000 LBS 




Figure 6-U: Symmetric Flutter Analysis - Mach 0.9 ~ Full Fuel Full Payload. 
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FREQUENCY, HERTZ FREQUENCY, HERTZ 


ANTISYMMETRIC FLUTTER ANALYSIS - CHORDWISE STIFFENED ARRANGEMENT 

MACH NO. = 0.9 
WEIGHT = 321,000 LBS. 




0 100 200 300 400 500 600 

VELOCITY, KEAS 


Figure 6-5: Antisymmetric Flutter Analysis - Mach 0.9 - 
Operating Weight Empty. 


ANTISYMMETRIC FLUTTER ANALYSIS - CHORDWISE STIFFENED ARRANGEMENT 

MACH NO. = 0.9 
WEIGHT = 750,000 LBS. 




VELOCITY, KEAS 


Figure 6-^6: Antisymmetric Flutter Analysis - Mach 0.9 - 
Full Fuel Full Payload. 
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FLUTTER SPEED , KEAS 


A'/ 

U 


20 MODES 


15 MODES 


10 MODES 




MACH 0.90 

WEIGHT = 321,000 LB. 
SYMMETRIC BOUNDARY 
CONDITION 


VELOCITY, KEAS 

Figure 6-Jt The Effects of Number of Modes on the Flutter Solution. 


SYMMETRIC BENDING AND TORSION MODE 


MACH 0.90 

WEIGHT ° 750,000 LB. 


NUMBER OF VIBRATION MODES USED IN THE FLUTTER ANALYSIS 

Figure 6-8: Flutter Speed Variation With Number of Modes 



probably could benefit from wing stiffening. However, the mode marked stabil- 
ity mode proved to be composed of the first elastic mode and rigid body modes 
and it was not expected to gain significantly in flutter speed by wing stif- 
fening. In view of its strong dependency on the rigid body stability char- 
acteristics, which during Task I were not well defined, this mode was not 
considered during the Task I optimization. 

Additional insight was gained by a flutter analysis of a configuration 
rigidized inboard from butt line 470 (Figure 6-9). The flutter mode with a 
flutter speed of 4l8 KEAS is very similar in character to the bending and 
torsion mode of Figure 6-4 with a flutter speed of 379 KEAS. 

On the basis of the results shown, and additional flutter results for 
other combinations of Mach number and weight , the symmetric bending and tor- 
sion flutter mode at M = 0.9 and 750,000 pounds was chosen for the optimiza- 
tion task. 


6.3 Optimization 


As indicated in the previous section the full fuel, full payload, 

750,000 pounds configuration, chordwise stiffened concept, at M = 0.9 was 
chosen for the optimization procedure. The flutter speed of the symmetric 
bending and torsion mode of the initial configuration was 379 KEAS. The flut- 
ter optimization task was to resize the structure to increase the flutter 
speed to 1-2 (468 KEAS) with a minimum weight penalty. 

To reduce the optimization task to a managable size relative to available 
resources the wing planform was .divided into eight regions (Figure 6-10). For 
each region a design variable related to the torsional stiffness and a design 
variable related to the bending stiffness were defined. Increasing the tor- 
sional stiffness was accomplished by increasing skin thickness and beam web 
thickness, which resulted in increased bending stiffness as well. 

Increasing the bending stiffness was accomplished by increasing the beam 
cap areas and resulted in a small increase in torsional stiffness due to dif- 
ferential bending. Design regions 9 and 10 are the bending stiffnesses of 
the engine supports. They were not varied during the Task I optimization 
procedure. 

The static reduction of the stiffness matrix from 1042 to l88 degrees of 
freedom led to non-linear relations between elements of the l88th order stiff- 
ness matrix and the design variables 3, as illustrated in a few examples in 
Figure 6-11. The stiffness matrix, as a function of all design variables 3, 
was approximated by: 

l6 

[K(3. )] * [K(0)] + y ([A.] + 3. [B ] + 3,? [C ]) (6.1) 

i=l 
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It was found that interaction between the design variables was minor so 


that the 


3. 3 . 


terms were neglected in the approximation. 


The optimization technique used is an interactive Computer Graphics pro- ^ 
gram, using the method of Incremented Flutter Analysis, in which successive 
coirfoinatiOns of design variables with decreasing total mass are generated 
while keeping the flutter speed exactly at the desired value. The effective- 
ness of design variables in increasing the flutter speed per unit mass is 
determined similar to the procedure described in Section 3.2. On the basis of 
the current values for the effectiveness the next resizing is determined by 
engineering judgment, after which new sensitivities are determined. The proc- 
ess is continued until the sensitivities have become approximately equal and 
no f\irther weight reduction is possible. If the number of design variables 
that play a role is small, say smaller than 20, this technique is very effi- 
cient in terms of resources required and elapsed time. 

The basic computer program in this optimization technique is executed 
with fixed modalization, although it can be interspersed at will with modal 
updating. The increase of the flutter speed to the desired value can also 
be done in several steps. The practical experience gained while following 
one particular sequence of steps is discussed. Reference is made to 
Table 6—2. 


The initial configuration (Sq) and its natural vibration modes obtained 
from a l88th order vibration analysis (VMq) had a flutter speed of 379 KEAS. 
The flutter equation was modalized with the lowest 20 of these modes. An 
optimization using fixed modes was performed and it was found that to reach a 
flutter speed of 468 KEAS, 1105 pounds per side had to be added. This first 
*’ optimized” solution (S^) consisted of 682.5 pounds of structure for the 
torsion design variable and 422.5 pounds for the bending design variable in 
Region 8 (Figure 6-10). For the configuration thus defined, a new vibration 
analysis (VMj_) using 188 degrees of freedom was conducted followed by a 
20 mode flutter analysis. The resultant flutter speed was 443 KEAS, 25 KEAS 
less than the desired value. Re-analyzing this structure using the original 
20 vibration modes led to the desired flutter speed of 468 KEAS, indicating 
that it was the use of fixed modes that caused the optimization process to 
result in a deficient structure. 

Using the newly computed vibration modes (VMp) in a fixed mode optimi- 
zation process indicated the need of an additional 492 pounds per side (for 
a total of 1597 pounds per side added) to reach a flutter speed of 468 KEAS. 
This second solution (S2) had l470 pounds of bending material and only 
126.7 pounds of torsion material. Thus, it appears that for twenty modes this 
process of fully optimizing and then performing a vibration analysis may not 
converge on the optimum distribution among the several design variables. This 
approach was considered unpromising and has not been pursued further. No 
updated flutter analysis of configuration S was performed. 
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Flutter Analysis 

Configu- 

ration 

Vibration 

Analysis 

188th 

Order 

Modes 

Used 

Order 

Flutter 

Speed 

KFAS 

Sq (base) 

"■'0 

™0 

20 

379 



VM. 

X 

20 

U43 



™0 

20 

i;6e 



VMq 

50 

UU5 

"3 

VM^ 

VI.l^ 

50 

U60 


Optimization 



Resultant Configuration 

Modes 

Used 



Weight Added to Base 
Configuration, lbs per Side 

Order 

f ication 

Region 

Torsion 

Bending 

Total 

c 

3; 

0 

^0 

Si 

8 

682.5 

U22.5 

1105 

VM^ 

20 


8 

— 

126.7 

1^70.0 

1596.7 


861.4 

82.0 

778.4 

1721.8 

704.0 

490.0 

755.0 

1949. 0 


Table 6-2; Arrow Wing Optimization 


























































To maximize the possibility of obtaining a meaningful optimized configu- 
ration the maximum number of modes,, 5Q, that could be accommodated on the 
interactive system was used. First configuration was analyzed, using 

50 of the VMq modes. This resulted in a flutter speed of ll45 KEAS, very 
close to the 443 KEAS found with 20 VM^ modes. Using the 50 VMq modes 
another optimization was performed, leading to S 3 with 86 l .4 pounds in tor- 
sion and“778.4 pounds in bending added in Region 8 and 82 pounds in bending 
in Region 7 . A l 88 th order vibration analysis of S3 resulted in the modes 
VM3, 50 of which led to a flutter speed of 460 KEAS. A second 50 fixed modes 
optimization was performed using the VM3 modes. This led to the configu- 
ration Si|: 704 pounds in torsion and 755 pounds in bending in Region 8 and 

490 pounds in bending in Region 7, for a total weight penalty of 1949 pounds 
per side. 

It was judged that remodalizing the Sl| configuration would lead to a 
50 modes flutter speed acceptably close to 468 KEAS and that the modes would 
change little, indicating that a minimum weight was approximated. 


6.4 General Conclusions 

Some of the practical experience with flutter optimization in -a realistic 
design environment gained during the Contract NASl-12288 flutter task is 
illustrated in the preceding section. Additional experience was gained dur- 
ing Task I- With two other structural design concepts, and during Task II with 
a mixed chordwise stiffened-monocoque design. Each optimization was done for 
a 750,000 pound configuration at M = 0.9. During Task II an optimization at 
M = 1.85 was added since the M = 0.9 optmization led to a flutter speed less 
than 1.2 Vj) at M = I .85 (Reference 19) . 

The practical experience gained by the present investigators during the 
entire arrow wing contract can be summarized in the following conclusions. 

1. An airplane with relatively complicated in-flight modal character- 
istics was optimized considering only one flutter speed at a time. 

2. The optimization of the arrow wing configurations was accomplished 
with the help of man-ln-the-loop techniques. The resulting' restric- 
tions regarding the number of design variables and structural and 
modal updating emphasize the need for more powerful optimization 
programs . 

3 . Reduction of the stiffness matrix and the associated non-linearity 
leads to the structural analysis being a dominant part of the cost 
of optimization. Future optimization studies should aim at reducing 
this cost. Such cost reduction may be found in using a simple struc- 
tural model or in using approximate methods. (See also discussion 

in Reference 1.) 


Frequent updating of the vibration modes assists in speeding up the 
attaining of a converged optimum design. An error in judgment, 
underestimating the frequency required, however, will be caught when 
the final, supposedly optimized design is remodalized and a check 
flutter analysis performed. Such a check should be a routine pro- 
cedure. The use of too low a number of vibration modes, or insuf- 
ficient design regions, however, is not "automatically" checked for 
accuracy, but would require effort outside the mainstream of the 
optimization procedure. 


T. DERIVATIVES WITH RESPECT TO DESIGN VARIABLES 

The resizing of the structure in a structural optimization with flutter 
constraints is directly related to the sensitivity of the flutter speed with 
respect to the design variables. In 1947, van der Vooren (Reference 20) pub- 
lished a method of expressing the flutter speed and mode of a slightly per- 
turbed configuration in terms of the flutter speed and mode of the original 
configuration and the perturbation. More recently, and directly in connection 
with flutter optimization work, Rudisill and Bhatia (Reference 6) have 
derived expressions for the derivative of the flutter speed with respect to 
design variables. A compact description of the method is included for com- 
pleteness in Section T-l- 

Section 7.2 is based on work performed during this study. It presents a 
method of obtaining derivatives of the damping at a given speed with respect 
to design variables. 

For both types of derivatives, the derivatives of the mass, stiffness 
and aerodynamics matrices are needed. Present finite element structural mod- 
elling technique leads to satisfactox’y structural representation for aero- 
elastic purposes while maintaining linear relationships between the mass and 
stiffness matrices and the design variables. Thus the mass matrix [M(mq)] 

arbitrary values of mp, can be written 


n . 

y m. [AM.] (7.1) 

i i 

i=l 


n 

i ^i 

i=l 


[AK.] 


(7.2) 


and the stiffness matrix [K(mi)], for 
as : 


[M(m^)l = [Mq] + 


and 


[K(m^)l = [Kq] + 
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where tK^] are base mass and stiffness matrices, respectively, 

and [AM^] and [AK^] are incremental matrices corresponding to a unit of 

m. i 

1 

If the order of the stiffness matrix is too large for efficient repeated 
vibration analyses, as required in an optimization procedure, static reduc- 
tion may be necessary; this may destroy the linear relationship shown in 
equation 7-2. However, as described in Reference 1, the derivative of the 
stiffness matrix for m^ = 0 can be obtained from an expression identical 

to equation 7.2. For each design variable for which m^ 7 ^ 0, new matrices 

[Kq] and [AK^] must be computed to determine the derivatives . 

The derivative of the aerodynamics matrix can be obtained by finite 
difference techniques or by the methods discussed in Section 5 of Reference 1. 


7.1 Derivatives of the Flutter Speed 


The flutter speed, V, and the associated reduced frequency, k = , 

are obtained from the characteristic flutter equation: 


_L [M] [K] - i p [A(ik)] 

c V 


{q} = 0 


(7.3) 


where g is the structural damping 


Since at flutter Y“0, equation 7-3 is valid for the k-method approach 
as well as for the p-k method approach. 

The matrix {q} is the modal column associated with the flutter solu- 
tion. For the same determinant, but associated with the transposed matrix 
problem, there exists the corresponding column {r} , Thus; 


L^J 



[M] k^ + [K] 

V 


•| p [A(ik)] 


0 (7.4) 


io6 


IS : 


The derivative of equation 7-3 with respect to the design variable m 


r 3M~] 2 

2 hm " 2 

P U -I P 


9k 2(l+ig) ^ p. 




9K 


9m 


(T.5) 


1 F-saI 

2 ^ [9kJ 


9m. 


{q} + 


- -i [M] k^ + [K] - I p [A(ik)] 

c V 


M- 


Pre-mult iplying 7.5 with L^J and invoking equation 7.U leads to: 
1 / 9M \ 


,2 2 ^ \ , 9k 2(l+ig) / -rr \ 

K - _ (rMq) ^ (rKq) 

C V 


(7.6) 


1+ig / 9K \ 1 / 9A \ 9k 


Equation 7-6 is a complex scalar equation. The scalar i* q stands 


9m 


for [rj {q}. Similarly the other triple products between parentheses 

stem from triple matrix products. In equation 7-6 k, V, {q}, and Li’J 
known solutions of the equations 7.3 and 7.4. Equation 7.6 represents two 


9k 


9V 


linear equations from which the two unknowns and “ can readily be 

obtained,. 

It may be of interest to return to equation 7.5 and determine the deriv- 
ative of the flutter mode: 




With and determined from equation 7.6 

9m 9m 


9V 


•M 


is the only 


unknown in equation 7.5. Since equation 7.3 defines only a distribution for 
{ 4)5 {q.) normalized. Let {q} be normalized to the last element, 

which then, by definition, is a constant. It follows that the last element 

of M is zero, which can be omitted if the last columns of the matrices 

pre-mult iplying m are omitted. Equation 7.5 then has one more equation 
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than there are unknowns and, arbitrarily, the equation corresponding 

to the last row of all matrices is omitted. ■ 

If [m|] means the last column of [M] is omitted, and [M] means 

the last row is omitted, the solution for can be written as: 

l.9mj 

T -1 


fe}-- 


- 4 ] -Ip [ A(ik) | 1 


V 


1 faM" 

2 |_ 9m 


k^ + 


_ ^fMlk ^ ^ 

9ra .^3 
c V 


3m 


r 1+ig [dK] 1 ^ f3Al 

^ LmJ ■ 2 " LisJ 


{q} 


(T.T) 


The complete column is given by: 




(7.8) 


The second derivative of the flutter speed is used in some methods of 
optimization that use a defined step size (Reference 1, Section 7). Two 


equations with the two unknowns 


9^k 


and 




can be derived by assum- 


9m. 9m. “““ 9m. 9m. 

1 J 1 J 

ing the derivatives in equation 7.5 are with respect to m. and then differ- 


entiating with respect to m. . Because of equation 7.1, 
. J 


2 

9 M 


9m. 3m. 
1 J 


= 0 . 


It can be shown that for T" 0 , 


2 

9 K 


9m . 9m , 
1 J 


= 0 . The unknowns 


9\ 


9m . 9m . 
1 9 


and 


9^V 7 

' can be expressed in terms of the originally known quantities and 

4..^ . , 4. ^ 3k 9k 3V 9V 

the quantities to be computed: ^ — » t — , , , . 

^ 9m. ’ 9m. ’ 9m. ’ 3m. ’ 9m. 

1 9 1 9 11 


3q 


and 


9q 


9m . 
9 
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7.2 


Derivatives of the Flutter Foot 

The flutter root, p, defines the frequency and damping of an in-flight 
mode at a given speed. It satisfies the characteristic equation of the 
flutter equation: 


-ViM] + (1+lg) [K] - i pv"" [A(p)] 


(q) = 0 


(T-9) 


Let p=(Y+l)h he a root of equation 7*9, and {q} the associated modal 

9y 

coliimn. The derivatives would he required in an optimization process in 

which there is a damping constraint. Damping may he a flutter constraint (as 
discussed in Reference 1, Section 7.3), hut could he a constraint related to 
airplane stability characteristics. 

can he obtained from = a+ih: 

3ra 9m 


9y _ a - 

9m k 


(7.10) 


^ is obtained from equation 7.9 as follows: 

8m 

2 2 12 

For a given Mach number and altitude V /c and — pV are constant, 

2 2 
V ~ ~ 

To simplify the notation — ^ [M] is replaced by [M] , (l+ig) [k] by [K] 

12 ~ ^ 

and pV"* [A(p)J by [A(p)J. Equation 7.9 then becomes; 


[M] p" + [K] - [A(p) 1 {q} = 0 


(7.11) 


Differentiation of equation 7.H with respect to the design variable 
m leads to: 


r 9M1 2 . r-. 

— p +2 LM] p ^ + 
L9fflJ 9m 


[I] - [I] 


9m 


{ql + 


(7.12) 


+ [M] p^ + [k] 
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Pre-multiplication by L^-J , vhere Ld satisfies equation T-13 for the 

3p . . 

flutter root of interest, leads to the expression for “ given in 
equation 


W[ 


[M] p^ + [K] - [A(p)] 


= 0 


(7-13) 


L^J 


3P _ 

3m 


M + p2 m 

|_3raJ ^ [_3mJ 


{1} 


(T.lH) 


L'J 


ra- 


2p [M] 


{q} 


Following a reasoning similar to that in Section T-lj the derivative of 
the flutter mode can he derived: 


-1-1 r 


{||}= - II + II - AIeII H + 2 [m] p 


3m 


(T.15) 


r?!! _ ( 

■■ 3 AI 9p 

1 

s 
1 

3pJ 3m 


{q} 


Assuming the derivatives in equation 7*12 are with respect to m^ , this 
equation can he differentiated with respect to mj and an expression 


for 


3% 


3m. 3m. 
1 J 


can he obtained, from which 


r.2 , a.h. - 2Yh.h. + a.h. 

9 y ^ c-yq _ 1 .1 1 j ,1 1 

3m. 3m. k ” ,2 

1 , J k 


(7.16) 


The a, h, c and d quantities in equation 7 -16 are defined by: 

S2p 


— = a + lb . 
3m. 1 1 

1 


= a . + ib . 
3m^ J J 


3m. 3m. 
1 J 


= c + id 


(7.17) 
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Since it is not certain whether acceptable approximations to [A(p)] can 
be derived from [A(ik)j given at discrete values , an alternate 

3y • • 

approach to determining is presented for consideration. 

In the p-k method (Reference 3) equation T-H is written as; 

|^[M] p^ + [K] - [A(lk)]j {q} = 0 (T.18) 

D,if f erentiation gives : 


■3M“| 2 ^ r-T 8p ^ [dK] FSAI 8k 

_taj P ^ '“IP to * [taj - [aid to 


{q} + 


(T.I9; 


+ [tH] + [ic] - [A(ik)] IIIIH 0 


]& 


M = ilk + (Y+i) ■' 

8m 8m ^ ^ 8m 


8k 


(7.20) 


Mult iplying'"^'q4at ion 7*19 ty I r J and substituting equation 7.20 leads 

87 ^ 8k 

to the following complex scalar equation with the unknowns and : 


/ 8M \ 

I" is V 


q\ p^ + (rMq)2pk 4^ + (rMq) 2p (y+i) 

oIU 


8k. 

/ 8K \ 


/ 8A \ 

' + 

8m 



V ^ V 


(7.21: 


As before, r q. is a short notation for the scalar quantity 

oUl 


I I 

L^J [_9m 


{q} , etc. 


No numerical evaluations have been made of methods to compute 


8m 


7.3 Derivatives of the Vibration Solution 

It is customary, and for practical reasons often mandaiory, to reduce 
the number of degrees of freedom of the flutter equation by modalization. In 
general the modalization is accomplished with the help of “vibration modes . 


Ill 


During an optimization process, in which the structure is resized in several 
steps, updating of the vibration modes is often necessary. Every update 
requires the solution of the vibration equation. 


Possibly the number of updates can be reduced if the change of vibration 
modes due to the changes in design variables is used. E.g., for the 
mode : 



(T.22) 


The procedure of Section 7-2 is followed. 


The vibration equation is; 


[- 


[M] 05^ + [K] {z} = 0 


(7.23) 


Eor a root m = o), : 
k 


[- 


[M] 0)^ + [K]j {z^} = 0 


and, because [M] and [K] are symmetric matrices: 


[- 


[M] + [K] = 0 


'] = 


(7.21+) 


(7.25) 


Differentiation of equation 7.24 with respect to m^ gives; 


-ft] 


2 *“k 

M, - tM] 


ocu r- 1 1“ 

^ {Z }-H - 

3m. 3m. k"^ 

1 L X-* L 


{z^} + - [M] + [K] 




(7.26) 


Premultiplication by ^z^^ leads to a solution for 


3u), 
k 

3m. 

1 


3ui, 


k 


am. P^JTmTR^ 


W f3M T 2 f3K ■ 

■ [am. J \ [3m 


{z^} 

k 


(f.27)'' 
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vibratiii°moSS reasoning in Sections 7.1 and 7.2 the derivative of the 
viDratixon mode can te denivedt 


{5}- - [• ["] • [“]J‘ [■ fe] ^ - w S • g] 


from vhich 




(7.28) 



f9z ^ 


k 


3m 

3rn hi 

1 

r“ij 1 

J 


A ne« modalization matrix is formed by assembling columns 


(7.29) 


f k }■ 

L newj 


^ vibration modes to be used,Tnto one 

3 2 

At any step in the optimization process the quantities defining 

3m. 

according to equation T.27 are ayallable. Howerer, the inyerse in ^ 

ation rt\e-rar?olS:„:!“^^^ 


7.4 Flutter Speed Derivatives with Variabl 


e Modes 


When modalizing the flutter equation with natural vibration modes these 
modes can be considered as a new set of coordinates suitable for dSinInf 

configurations similar to the configuration 
speed and flutte determined. In determining the flutter 

tL formi.fSon of modes are then considered constant, 

as^umnt?^ ni ^ ^ J derivatives in Sections 7-1 and 7-2 is based on the 

assumption of constant modalization matrices. 

^ 4 - lioweyer, the use of updated vibration modes is considered an inteeral 

stant° modalization matrix cannot be considered con- 

be inclSdefl^ AerlTatlves of the modalization matrix must 

resSti ofkctL" in Reprance SI. A procedure, based on Se 

7-3j IS outlined in the following paragraphs. 
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r9M“i . 

[arj 


Consider I :: — I in equation 7- 5- If equation 7*3 is a modali zed flutter 
equation, [M] is a modalized mass matrix. 

[M] = [M] [7] (7.30) 

where [m] is the basic matrix and [z] is a matrix. modal columns. 

It follows: 


raMl r-T-i raM“| 


[z] + 


_T 

9z 

am 


■I 1 mi 


[M] [z] + [z ] [m] 


az 

am 


(7.3i: 


The second and third term on the right hand side of equation 7.31 are 
extra terms due uo letting the modes be functions of the design variables. 

An equation similar to equation 7-31 can be written for 1—1 

The expression for I I becomes: 


9m 


aA p — T 1 3A”1 

L^J - ^ [stj 


tzl II . [A] 



3z‘ 


am 


[A] [7] 


(7.32) 


The matrix 


8z 

IX •— 

dm 


is formed by assembling columns 


1 3m/’ 


given by 


equation 7 . 29 , for as many |z I as there are vibration modes used in the 
modalization. i 


Inclusion of 



in the formulation of the flutter speed and flutter 


root derivatives is expected to lead to more accurate derivatives and con- 
sequently to fewer resizing steps in the optimization process. This, 
presumably, will reduce the number of vibration analyses required for remodal- 
ization, which is the same result as is expected from the procedure presented 
in Section 7.3. 

That a potential for an increase in efficiency due to the inclusion of 

razl . 

3m I exists, follows from the following consideration. 

In choosing design variables, practical considerations may favor local 
stiffnesses as independent design variables. To reflect the effect of a local 
stiffness variation in the modalized flutter equation, the ^’resolution power” 
of the modes must be sufficient to recognize this stiffness variation. When 
the flutter speed and flutter 'root derivatives are determined under the 


llH 


assumption of constant modes it seems self-evident that more modes are 
required for derivatives of a desired accuracy than when the modes are con- 
sidered functions of the design variables. In the latter case, if the 
resolution of the unmodalized vibration equation is sufficient to recognize 
local stiffness variations, then this resolution power' is partially 


preserved through the inclusion of the 


terms. 


Numerical investigations to confirm the above expectations are 
recommended.- 
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